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1 Introduction 


This report documents the first stage of research into Control Law — Control 
Allocation Interactions. A three-year research effort was originally proposed: 

1. Create a desktop flight simulation environment under which experiments 
related to the open questions may be conducted. 

2. Conduct research to determine which aspects of control allocation have 
impact upon control law design that merits further research. 

3. Conduct research into those aspects of control allocation identified above, 
and their impacts upon control law design. 

Simulation code was written utilizing the F/A-18 airframe in the power- 
appraoch (PA) configuration. A dynamic inversion control law was implemented 
and used to drive a state-of-the-art control alloction subroutine. 


2 Simulation 

The airframe used was derived from the F/A-18 model already implemented 
in CASTLE. The airframe is not realistic, but is intended to be a test-bed for 
further research. The greatest area in which the test-bed simulation differs 
from the original airframe is in the treatment of control deflections. There are 
essentially two sets of control effectors: 

1. The original control effectors in the F/A-18 airframe model. These are 
used only for initial trim and subsequent scheduling. 
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2. A duplicate set of control effectors that have linear effectiveness. This 
control set is the input to the control allocator, and the forces and moments 
they generate are superimposed on those of the bare airframe and original 
control set. Rate limits of the duplicate set are the no-load rate limits of 
the original controls. Position limits of the duplicate control effectors are 
referenced from the trim or scheduled postions of their counterparts in the 
original controls. 

The rationale for incorporating a duplicate control set was to provide a 
constant, linear control effective matrix with flexibility for future variations and 
modifications. The control deflections are calculated for the trimmed/scheduled 
flight condition using the original F/A-18 nonlinear table lookups. The control 
deflections calculated from the allocator to produce desired moments use the 
control effectiveness matrix obtained from linearizing the F/A-18 aero database. 

2.1 Simulation Code 

There are six files that are used in the simulation of the airframe: Aero.f, 
Aeropa.f, Control. f, Constants. f, Engine. f, and Alloc. f. 

2.1.1 Aero.f 

The aero code first calls Aeropa.f to calculate the aerodynamics of the sched- 
uled/trimmed flight condition. The code then combines the aerodynamics from 
the non-linear scheduled/trimmed flight condition and from the control deflec- 
tions calculated in the allocator to produce the desired moments. 
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2.1.2 Aeropa.f 


This code is taken from the F/A-18 simulation and modified slightly to include 
the control positions that are used in the table lookups later in the code. This 
is the only code that gives the airframe F/A-18 like characteristics. All added 
code is at the top of the RUN section. 

2.1.3 Control. f 

Stick and rudder pedal commands are taken as inputs and converted into an al- 
pha command a cm d, beta command 0 cm d, and a roll rate command p cm d These 
commands are input to a simple dynamic inversion control law that generates 
desired moments for the control allocation subroutine. First, a cm d and 0 crn d 
are converted to desired accelerations Wdes and Vdes- 

U'cmd ~ tan Qcmd 
l^des — ^cmd) 

Ucrnd ~ ^ $in 0cmd 
Udes = ^ v Vcmd) 

Next, ubdes and Vdes are applied to inversions of the body-axis force equations 
(treating q and r as controls): 

tides +P L ’ — g cos 9 cos <£> — Z/m 

Qcmd — u 

— I ^de.s+pu> + '? cos 0 sin <p+Y / m 

r cmd — u 

These two inversions are made as perfect as possible by using actual aircraft 
states, and the last calculated values of the body-axis forces Y and Z from the 
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aerodynamic calculations. First-order responses are specified for the desired 
angular accelerations, 

Pdes — ^ p ( P ~ Pcmd ) 

Qdes = ^ q {Q Qcmd ) 

^ des — ^ r ^cmd) 

Finally, the desired body-axis moments, required to obtain the desired ac- 
celerations, are calculated from inversions of the body-axis moment equations: 

/'"'♦C i Ixxpden “ ^xz ^de9 + ( ?z z ~ ly y )Q r ~ Ixz PR 

+ qSb 

nc _ r^a , tyy<i<ies+(Ij:x-Izz)p T +Ixz(p 2 ~r 2 ) 

+ ^Sc 

C f~^ a ~ jjg * pdf ? [ z z 1~de s )P < ?~^~ £l * Q r 

^n des — L 'n f qSb 

The moment coefficients C“, C^, and are the last calculated values of 
the body-axis moments. Since control-generated moments are superimposed 
on these values, they are the moments generated by the bare-airframe plus 
scheduled control deflections. The trimmed flight control deflections are used 
to calculate moments for the current flight condition to be used in the restoring 
algorithm. The attained moments are calculated next using the control deflec- 
tions from the last time step for comparison purposes with the desired moments. 
The desired moments, along with the required inputs, are input to the allocator 
to produce the required control deflections. The last step is to check the control 
deflections against the limits and reset them accordingly. 

2.1.4 Constants. f 

This section of code sets the model specific contants. 
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2.1.5 Engine. f 


The engine model is taken from the Stevens &: Lewis F-16 model [1]. 

2.1.6 Alloc. f 

This code is the control allocator that produces required control deflections for 
desired moments. This code is explained in detail later. 

3 Desktop Simulation 

The F/A-18 PA model was first implemented on the UNIX-based CASTLE. 
The conversion of the simulation to the desktop PC required the CASTLE 
offline help menu provided with the PC version of CASTLE. Some additional 
steps were taken to complete the compilation of the airframe. The steps are as 
follows: 

1. The directory structure from UNIX was copied to the CASTLE airframes 
folder. 

2. A project was created in Microsoft Studio 6.0 foilwing the F/A-18 project 
already included with the PC version of CASTLE. All custom builds were 
set up in the same way the F/A-18 project had them set up 1 . The custom 
builds were implemented on symbols. sdf and all the FTP data files. 

3. In symbols. sdf the realtime CDF section was changed to resemble the 

l The offline CASTLE help explains a different way of setting up custom builds, but did 
not work. 
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F/A-18 realtime CDF section in the corresponding symbols. sdf. The 
reason is just a difference in structure between the PC CASTLE and 
UNIX CASTLE. 

4 Control Allocation Algorithm 

4.1 Introduction 

The control allocation algorithm is a FORTRAN implementation of the bi- 
secting, edge-searching algorithm. The theory behind the allocation code is 
explained in detail in [2]. Following is a step-by-step explanation of the code. 
Line numbers correspond to those in the attached file “Alloc. f’. 

4.2 Subroutine DA3 

4.2.1 Diagnostics 

The sections of code that depend on the DIAGS flag are debugging tools that 
can be used to dump several relevant variables. Because a great volume of 
output is generated, the DIAGS flag should be used sparingly. 

4.2.2 Code Description 

Lines 0126-0146: Array CSPHI is a table of sines and cosines of angles, 
beginning at 45° and proceeding through 20 bisections. 

Lines 0191-0208: The desired moments are checked for zero length, 
and a vector of zero control deflections returned if they are; otherwise the 
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vector is normalized. 


Lines 0210-0219: The initial rotation is performed using the transfor- 
mation generated by subroutine DCGEN to align the desired moment 
direction y 3 d with the yi axis. Subroutine DCGEN is an implementation 
of the method described in [2, Section 5.1]. Lines 0212-0219 perform the 
matrix multiplication, B$ = T£? 3 orig . 

Lines 0221-0231: The controls that generate the moment with the max- 
imum yi component are found by examining the sign of the first row of B 
and setting the control to its maximum or minimum, depending on that 
sign. The controls are first set to ±1 (object notation) and then set to 
their actual limits by subroutine SETU. 

Lines 0233-0243: This section of code was added to deal with the finite 
precision of computer math. The variable TOL is a distance in moment 
space that is related to the smallest bisection angle to be used, at the 
distance from the origin of the vertex just determined (maximum y\ com- 
ponent). TOL is used in subsequent code to resolve near-zero numbers. 

Lines 0264-0265: Subroutine R20 solves the 2-D problem for the pro- 
jection onto the current y\-y 2 plane. R20 returns the object-notation 
vector of controls of the intersecting edge in variable Ut, the control that 
defines that edge in variable IU , and a ±1 value in variable INFRONT 
that is -hi if the edge is in front of, and —1 if it is behind the yi- 1/2 plane. 
The three variables TEMP2, TEMP3, and TEMP4 contain respectively 
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the sorted list of controls (ITHETA) with an additional zero between the 
two controls at the ends of the intersecting edge, the number of vertices 
in the list (NANGS), and the index of the position in the list of the addi- 
tional zero (INDX). Finally, the logical variable IS VERTEX signals that 
the desired moment points directly at a vertex. 

Lines 0274-0293: This section of code has no counterpart in reference 
[2]. It was added during debugging and found to improve the success rate 
of the algorithm (decrease the number of estimations required). The most 
recently found edges that were in front (Last In Front, LIF) and behind 
(Last In Back, LIB) are saved. Theoretically the last two edges found will 
be LIF and LIB, but in some cases they were not. 

Lines 0295-0299: If R20 reports a vertex in ISVERTEX, the controls 
that determine that vertex, and the saturation of the desired moment, are 
calculated by a call to DO VERT, and the subroutine is exited. 

Lines 0304-0322: This section of code initializes several variables, in- 
cluding the rotation matrix T22. 

Lines 0333-0510: This is the main loop, in which the 2-D problem is 
repeatedly solved for different rotations about the y x axis. 

Lines 0335-0342: Used during debugging, retained for possible 
future use. 

Lines 0344-0349: Rotation about y x . B 1 is the operative B matrix 
throughout. Code performs operation T ■ B. 


9 



Lines 0360-0364: The last returned values of ITHETA, NANGS, 
and INDX are assigned to those variables to be saved when TEMP2, 
TEMP3, and TEMP4 are overwritten by R20. 

Lines 0366—0367: Call to R20 to solve the 2-D problem for the 
current orientation of B1 about the t/i-axis. 

Lines 0376—0395: The edge identified by R20 is assigned to LIF 
or LIB according to the sign of the variable INFRONT. 

Lines 0397-0401: Another vertex check. 

Lines 0411-0495: Executed when the most recent and the prior 
edges differ in sign of their 1/3 component, as indicated by the vari- 
ables INFRONT and WASINF. This section of code is the implemen- 
tation of the description given in [2, Section 5.3]. Through line 0436 
the code is doing housekeeping and (possibly) diagnostics. 

Lines 0438-0457: This section reflects a subtlety in the imple- 
mentation of the algorithm not described in [2]. The prior edge 
was identified using a different B matrix than the most recent 
edge. All relevant information regarding the prior edge is con- 
tained in the saved variables ITHETA, NANGS, and INDX. At 
lines 0456-0457 a call is made to subroutine GETEDGE, which 
is also called as the last step of subroutine R20. 

Lines 0459-0478: More last-in-front, last-in~back checking, 
and lines 0480-0484 deal with vertex checking. 
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Lines 0486-0493: Check the last two edges identified to see if 
they comprise the solution facet. If they do not, the LIF and 
LIB edges are checked. Both sets of edges are checked using sub- 
routine ISFACET, described below. Output from ISFACET 
consists of the logical ISOK, numbers of the two defining con- 
trols in IUOUT and JUOUT, and controls (in object notation) 
at three vertices of the facet as columns of the array U123. 

Lines 0518-0519: If the variable ISOK is false, the correct facet has 
not been determined and the maximum number of bisections has been 
performed. One last check of LIF and LIB is performed. 

Lines 0520-0574: If ISOK is true, the solution is calculated. Otherwise 
(lines 0572-0573) the solution is estimated. 

Lines 0521-0565: A straightforward implementation of [2, Equa- 
tions (13) and (14)]. M123 is the matrix [e 3( i (v\ - v£) (Vj - V3)] 
in [2, Equation (13)]; variables AA, BB, and CC correspond to c* 3 , 
C\ t 2, and Ci <3 respectively; and MTEMP is v^. The variable UDA is 
the same as u u in [2, Equation (14)], except that it has been scaled 
as necessary. 

Lines 0572-0573: The estimator is called. 
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4.3 Subroutine DCGEN 


This subroutine is a straightforward implementation of the initial transformation 
algorithm described in [2, Section 5.1]. 

Lines 0789-0798: The desired moments are normalized using double 
precision math. 

Lines 0810-0824: If one or more of the leading components of the nor- 
malized moment vector are zero, the size of the problem is reduced. 

Lines 0829-0833: The first row of the transformation matrix is set to 
the normalized desired moments. 

Lines 0837-0850: The remaining terms are calculated in the three 
nested do-loops in [2, Equation (4)]. 

Lines 0858-0868: The last section of DCGEN ensures that the deter- 
minant of the transformation matrix is -hi. 

4.4 Subroutine R20 

To find the edge that the desired moments direction is pointing to, the subrou- 
tine R20 is implemented. The theory behind this subroutine is in [2, Section 
5.2.2], All calculations in this subroutine are done in the yi-y 2 plane. 

Lines 0928—0961: The y 2 component of the point with the maximum 
t/i component (UMAX in object notation, XUMAX in control notation) 
is calculated to determine its sign. The desired moment is checked to see 
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if its direction points towards a vertex of the attainable moment subset. 
If it is a vertex the subroutine is exited and the allocation carries on. 

Lines 0976-0992: Implementation of [2, Step 1, Page 20]. The array 
THETA is the needed part of the set £$, and ITHETA that of £ u . Once 
the angle is found, rr is added or subtracted from it if the absolute value 
is greater than tt/2 and depending on the sign of the angle. In this way, 
the angles of just the vertices with positive y\ components are generated. 

Lines 0998-1010: The angles are sorted in a clockwise or counter- 
clockwise manner starting with the vertex that has the largest \j 2 com- 
ponent. The manner in which they are sorted depends on the sign of the 
y 2 component of the maximum vertex, recorded in SY. 

Lines 1013-1025: A zero is inserted in THETA and ITHETA to mark 
the point at which the angle changes sign. 

Lines 1034-1036: THETA, ITHETA, and NANGS (the number of an- 
gles generated) are sent to subroutine GETEDGE to finish the solution 
to the 2-D problem. Subroutine GETEDGE is provided separately so 
that it could be called independently from DA3, as described above. 

4.5 Subroutine GETEDGE 

This subroutine is part of the explanation in [2, Section 5.2.2]. 

Lines 1090—1127: The first loop in this subroutine is looking for a sign 
change in the y 2 component between ordered vertices. Since the vertices 
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were sorted in the manner described, the solution edge will be the first 
one encountered in traversing the edges starting with the first vertex. The 
list is stepped through in the proper direction by the index IX = IX-SY. 

The previous ij 2 value is stored before the next y2 value is calculated. This 
new value is compared to the previous one determining whether the edge 
crosses the y\ axis. If the do loop continues, U2 is set to the next vertex 
by changing the sign of the control that is defining the current edge. The 
index is updated accordingly with the sign of y? and the process starts 
again until the edge is found. The do loop is exited when a new point is 
found that has a different sign than the point before. 

Lines 1129-1181: This section deals with possible failure of the previous 
loop to find an edge, as indicated by (SY.EQ.SSY). The starting values 
of relevant variables are restored, and the vertex list is traversed in the 
opposite direction. The first loop should always find the proper edge when 
GETEDGE is called from R20, but the first loop may fail when called 
from within DA3. The list is traversed in the opposite direction by the 
index IX = IX+SY. Implementation of this section of code was the reason 
for inserting a zero in the ordered list of vertices. 

Lines 1190-1200: One or the other of the previous two loops will have 
identified U2 (a vertex in object notation) and JU (the number of the 
control that defines the edge). U2 is converted to control notation using 
the subroutine SETU. The third row of the B matrix is applied to the two 
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vertices that define the solution edge to determine the component in 
moment space at the point where the t /2 component of the edge is zero. If 
the y 3 component is positive, the edge is described as “in front”, whereas 
if the y 3 component is negative, the edge is ‘‘behind” the line defined by 
the direction of the desired moments ^ 3 . 

Lines 1202-1217: A final vertex check is made and the subroutine is 
exited. 

4.6 Subroutine DOVERT 

Lines 693-728: If it was determined that the desired moments points directly 
to a vertex the subroutine DOVERT is called. DOVERT uses the maximum 
or minimum controls that make up the vertex and calculates the total moment 
from there, scaling it appropriately. The allocator subroutine is then exited and 
the simulation carries on. This case is rare during simulation, but may occur. 

4.7 Subroutine EST 

The theory behind the estimator subroutine is explained in [2, Section 5.4.2). 
Lines 0604-629 The subroutine starts with the last tw*o edges that the allo- 
cator had found and creates a facet by setting the appropriate control to — I or 
+ 1. SETU is used to assign actual control limits to these points which are then 
put into moment space using the control effectiveness matrix. 

Lines 631-669 An interpolation is then made with the estimated facet 
vertices to determine the solution. 
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Lines 671-686 The moments are calculated using the estimated control 
positions and then scaled with the saturation limits. 

4.8 Subroutine ISFACET 

The subroutine is used to test the facet found by DA3. The subroutine uses 
the two defining controls from DA3 to find a facet from scratch that these two 
controls define. This algorithm is the subject of reference [4]. 

Lines 1236-1251: Zeros are set in the appropriate positions of the vertex 
arrays so that two edges are defined for the facet. The dimension of the 
union (see (2, Section 4.2]) of the two edges is determined. If the union is 
not two dimensional, then the edges can not form a facet; ISOK is set to 
false and the subroutine is exited. 

Lines 1253-1320: For the two dimensional case the routine begins to 
calculate from scratch the facet that is determined by the two defining con- 
trols. The method used is completely independent of the edge-searching 
method and is explained in [4]. 

Lines 1255—1287: This section of code was lifted from earlier FOR- 
TRAN implementations of the facet-searching allocation method de- 
scribed in reference [4]. The facet defined by the two controls is in 
the variable TESTFACET. 

Lines 1291-1311: The facet TESTFACET is compared with the 
object OBJ that was generated by R20. If they are different, the 
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facet opposite TESTFACET (also generated by the same two con- 
trols) is tested (lines 1300-1311). 

Lines 1322-1336: If the facet just found is the same facet as the one 
that was found from the allocator, then U123, which is the matrix whose 
columns correspond to controls that generate three of the verticies that 
make up the solution facet, is assembled and returned. 

4.9 Miscellaneous Subroutines 

4.9.1 MINNORM 

The purpose of the minimum-norm restoring solution is to keep the controls as 
close to their trimmed control position as possible. The usual minimum-norm 
solution keeps the controls as close to zero as possible, however, in this applica- 
tion the zero position is redefined as the trimmed/scheduled control positions. 

Lines 1496-1531: The subroutine is started by finding the total control 
position for the current time step and calculating the total attained moment. 
Lines 1533—1554: If the control limits are zero, the routine is returned and 
no restoring takes place. Otherwise, the difference between the pseudo- inverse 
solution redefined at the trim condition, and the controls given by the allocation 
routine are used to find a delta control position that will drive the controls 
towards the trimmed position. This delta control position is scaled according 
to the control limits and a new restored control position is returned. 

For more information on control restoring, refer to Bolling. [3, Ch. 4] 
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4.9.2 SORTC 


Lines 1351-1492: A sorting subroutine downloaded from the National Insti- 
tute of Standards and Technology (NIST) GAMS (Guide to Available Math- 
ematical Software) at http://gams.nist.gov/. This particular algorithm was 
chosen for its efficiency, and for the fact that it returns a sorted index vector 
along with the sorted vector. 

4.9.3 INVMAT3 

Lines 0740-767: A brute force matrix inversion subroutine. Good only for 
3x3 matrices. 


5 Verification Data 

Sample runs are included to verify the airframe. The four tests cases used are a 
trimmed flight condition, a step in the longitudinal stick, and step doublets in 
the lateral stick and rudder pedals. The MANGEN command in CASTLE was 
implemented to produce the desired stick commands. Complete MATLAB files 
of the four cases are attached as trinudecl I. mat, long.decll .mat, lat.decll.mat, 
and dir.decll .mat. 

The plots include selected states of the airframe along with the trimmed/scheduled 
control positions and the allocated control positions. 

Figure 1 shows the time histories of the six global controls in a trimmed flight 
condition at 8.1 degrees angle of attack, 1200 ft, and 231.52 These settings 
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are the default when the airframe is loaded. Some settling of the controls to 
achieve steady state is noted. 

Figure 2 shows time histories for a step input in longitudinal stick of 2.5 
inches aft from center. The airframe was initialized to the trim conditions 
described above and the stick step implemented at time = 1 sec for 1 second. 

Figure 3 shows the time histories for a step douplet in lateral stick. The 
lateral stick was driven right 2 inches from center at time = 1 sec for 1 second 
and then left of center 2 inches for 1 second. 

Figure 4 shows the time histories for a step douplet in rudder pedals. The 
pedals were driven right 2 inches from center at time — 1 sec for 1 second and 
then left of center 2 inches for 1 second. 
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0058 C 

0059 C NONE 

0060 C 

0061 

0062 SUBROUTINE DA3 (UDA, SAT, IERR, 

0063 C AS FUNCTIONS OF 

0064 5c B, MDES , U_MIN, U_MAX, M, NBI , TIME, DIAGS) 

0065 

0066 C 

0067 C 

0068 C DECLARATION SECTION 

0069 C 

0070 C 

0071 IMPLICIT NONE 

0072 C ** Parameters 

0073 

0074 

0075 C ** INPUTS: 

0076 

0077 INTEGERS IMODE 

0078 

0079 C ** OTHER LOCALS: 

0080 

0081 BYTE CONPAR, CTLBUF 

0082 LOG I CAL* 4 DIAGS, DIDSWITCH, INITIALIZED, ISOK, ISVERTEX, STUCK 

0083 INTEGERM I, ICOUNT, IERR, INDX, ITHETA (21 ) , IU, IUOUT 

0084 INTEGERM IUTEMP{20), I_LIB, I_LIF, J, JU, JUOUT, K, M 

0085 INTEGERM MAXSTEPS, NANGS , NBI, NMAX , STEPS, SY, TEMP2(21) 

0086 INTEGERM TEMP3 , TEMP4 , Ul(20), U123 (20,3), U2(20), UMAX ( 2 0 ) 

0087 INTEGERM U_LIB(20), U_LIF{20) 

0088 REAL *4 MINV{3,3), MTEMP ( 3 ) , MXMAX ( 3 ) , COSPHI, CSPHI(2,20) 

0089 REAL *4 ABC(3), NORM, PI, SAT, SINPHI, DET, AA, INFRONT, T{3,3) 

0090 REAL *4 T22(2,2), BB, BTEMP (2 ) , CC, TIME, TOL, TOLANG, B(3,20) 

0091 REAL *4 M123(3,3), MAXNORM, UDA(20), Bl(3,20), MD(3), MDES(3) 

0092 REAL *4 U_MAX(20), U_MIN(20), WASINF, XU123(20,3), XUMAX(20) 

0093 REAL * 4 XUTEMP(20), Y 
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0097 C COMMON SECTION 
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0100 
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0103 C 
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0115 C 

0116 C 

0117 C DATA SECTION 

0118 C 

0119 C 

0120 

0121 DATA INITIALIZED/. FALSE./ 

0122 DATA PI/3.141592653589793/ 

0123 C 

0124 C Table of cosines and sines of bisection angle 

0125 C 

0126 DATA CS PHI/ 

0127 & 7 . 07 1067 8 1 18 6 54 7 5e- 01, 7 . 0710678118 6 5476e - 01 , 

012 8 Sc 9. 23879532511286 7e -01, 3.826834323650898e-01, 

0129 & 9 . 807852 8 04 03 23 04e- 01, 1. 95090322016 12 82e-01, 

0130 Sc 9. 95 18 47266 72 196 9e -01, 9. 80 17 14 03295606 Oe- 02, 

0131 Sc 9 . 98795456205 1724e- 01, 4. 90676 743274 18 Ole- 02, 

013 2 Sc 9 . 9 96 98 8 18 6 962 Q42e- 01, 2 . 4 54 122 8522 9 122 9e - 02 , 

0133 Sc 9 . 99924 7 0183 9 14 4 5e- 01, 1.22 7153 82 8 571993e-02, 

0134 Sc 9 . 99981175282 6 Olle- 01, 6 . 13 5 88464 9154475e - 03 , 

013 5 Sc 9 . 999952938095762e-01, 3 . 0 67 9 56762 96 5 976e- 03 , 

013 6 Sc 9 . 9999 8 82 345 170 19e- 01, 1. 5339801862 84766e- 03, 

013 7 Sc 9 . 99999705862 8822e-01, 7 . 66 9 903 18 7427 045e- 04 , 

0138 Sc 9. 99999 92 646 57 17 9e -01, 3.834951875713956e-04, 

013 9 Sc 9. 99999981616429 3e -01, 1. 9 17475 973 10703 3e- 04 , 

014 0 Sc 9 . 999999954 0410 7 3e - 01, 9 . 58 73 7 9 90 95 97734e- 05 , 

0141 Sc 9. 999999 98 85102 6 9e -01, 4. 79368996030668 8e-05, 

0142 Sc 9 . 999999997127567e-01, 2. 396844980 84 1822e-05, 

014 3 Sc 9. 99999999 92 8189 2e -01, 1.1984224 90506971e-05, 

0144 Sc 9. 99999999982 0472e -01, 5. 9 92 1124 5264242 8e- 06 , 

014 5 Sc 9. 99999999995511 8e -01, 2. 9 96 05 62263 34 66 le- 06 , 

0146 Sc 9. 99999999998878 0e-01, 1.498 02 8113169011e-06/ 

0147 

0148 C 

0149 C 

0150 C INITIALIZATION SECTION 

0151 C 

0152 C 

0153 

0154 IF ( (IMODE.LE. -2) .OR. . NOT . INITIALIZED ) THEN 

0155 ENDIF 

0156 C- - 

0157 C 

0158 C RESET SECTION 

0159 C 

0160 C 

0161 

0162 IF ( (IMODE.LE. -1) .OR. ( .NOT. Initialized) ) THEN 

0163 Initialized = .TRUE. 

0164 C 

0165 C IERR = 0 FACET FOUND, ABC OK 

0166 C IERR = 1 FACET NOT FOUND, INTERPOLATED SOLUTION 

0167 C 

0168 ENDIF 

0169 C 

0170 C 

0171 C RUN SECTION 
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0172 

C 




0173 

C-- 




0174 

c 




0175 


IF (DIAGS) THEN 



0176 


WRITE (*, 1 i/ASC! ' ) ’ En 

cerir.g DA 

5 

0177 


WRITE (*,' ;A50? 1 ) 1 Cal 

1 ir.g argur 

ent s 1 

0178 


WRITE (*, • 'A: , (Z19 .6} 1 

) 1 *DA3 * 

3 { 1 , : J 

0179 


WRITE ( * , 1 ..A3C,6Sia .6) 1 

) ' * DA3 * 

3 ( 2 , : ) 

0180 


WRITE (*, • ;A3C,6E13.6) ' 

) ■ * DA 3 * 

3 ; 3 , : ) 

0181 


WRITE ( * , ' { A 3 0 , 3 218. 6) 1 

) ' * DA 3 * 

MDES 

0182 


WRITE ( * , 1 *: A 3 0 , 6 F ’ 4 . 6 ) 1 

) ’ *CA3 * 

IJ MIN 

0183 


WRITE ( * , 1 ( A3C , 6? 14 . 6 ) 1 

) 1 * DA3 * 

U MAX 

0184 


WRITE (*, ' ' A3C , 13) ' ) 

' * CA3 * 

M 

0185 


WRITE (*, ' A 1C , 13 ' ' ) 

1 * D + 

mb : 

0186 


WRITE (*, 1 • A3C , FI = .6) ' ) 

t * n u * r ’ ' r 

M E = 

0187 


ENDIF 



0188 





0189 


INFRONT =1.0 



0190 





0191 


NORM =0.0 



0192 


DO 1=1,3 



0193 


NORM = NORM + MDES ( I ) *MDES ( I ) 


0194 


ENDDO 



0195 

c 




0196 


IF (NORM .EQ. 0.0) THEN 



0197 


IERR = 0 



0198 


SAT =0.0 



0199 


DO I = 1 , M 



0200 


UDA ( I ) = 0.0 



0201 


ENDDO 



0202 


RETURN 



0203 


ENDIF 



0204 





0205 


NORM = SQRT (NORM) 



0206 


DO 1=1,3 



0207 


MD ( I ) = MDES (I) /NORM 



0208 


ENDDO 



0209 

c 




0210 


CALL DCGEN (T, MD) 



0211 

c 




0212 


DO I = 1,3 



0213 


DO J = 1,M 



0214 


B1 (I, J) = 0.0 



0215 


DO K = 1,3 



0216 


B1 ( I , J) = B1(I, J) + 

T (I , K) *B(K, J) 

0217 


ENDDO 



0218 


ENDDO 



0219 


ENDDO 



0220 

c 




0221 


DO I = 1,M 



0222 


IF (B1 (1, I) .EQ. 0 . ) THEN 



0223 


UMAX ( I ) = 0 



0224 


ELSEIF (Bl(l, I) . LT .0.0) 

THEN 


0225 


UMAX (I) = -1 



0226 


ELSE 



0227 


UMAX ( I ) = 1 



0228 


ENDIF 




(B ( 1 , I ) , 1=1, M) 

(B ( 2 , I ) , 1 = 1, M) 

', (B ( 3 , I ) , 1 = 1, M) 

', (MDES(I), 1=1,3) 
', (U_MIN ( I ) , 1=1, M) 

1 , ( U_MAX ( I ) , 1 = 1, M) 

' , M 
' , NBI 
TIME 
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0266 C 

0267 

0268 

0269 

0270 

0271 

0272 

0273 C 

0274 

0275 

0276 

0277 

0278 

0279 

0280 
0281 
0282 

0283 

0284 

0285 


Friday, March 23, 2001 12:46:16 PM 


Page 5 of 28 
Printed For: Bull Durham 


ENDDO 

CALL S ETU ( XUMAX , UMAX , U_M I N , U_MAX , M ) 

TOLANG = CSPHI ( 2 , MIN < 2 0 , 2 *NBI ) ) 

DO 1=1, 3 

MXMAX (I) =0 . 

DO J= 1 , M 

MXMAX (I) = MXMAX (I) + B 1 (I,J) * XUMAX ( J ) 
ENDDO 
ENDDO 

MAXNORM = SQRT (MXMAX (1) *MXMAX(1) 

Sc + MXMAX ( 2 ) * MXMAX ( 2 ) 

Sc +MXMAX ( 3 ) * MXMAX ( 3 ) ) 

TOL = MAXNORM* TOLANG 

IF (DIAGS) THEN 


WRITE ( * , 
WRITE (*, 

' /A5C } ' ) ' Frail 

' A 3 0 , £13 . S ) ' ) 

ml nary 
*DA3 * 

J a _ c s ' 

_ 


NORM 

WRITE (*, 

' A30, 3 F14 . 6) ' ) ' 

* DAS * 

yo 

= 


( MD ( I ) , 1=1,3) 

WRITE ( * , 

' A3C, 3?14 .6) ' ) ' 

*DA3* 

T '■ 1 : ! 

= 


(T ( 1 , I ) , 1=1,3) 

WRITE (*, 

’ ■ A3 0, 3 F14 . 6) ’ ) 1 

*DA3* 

T ( 2 , : ) 

= 


(T (2 , I) , 1=1,3) 

WRITE (*, 

• ;A30, 3F14.6) ' ) ' 

*DA3 * 

: ' 3 , : ) 

= 


(T(3 , I) , 1=1,3) 

WRITE (*, 

' : A3 0 , 6E13 . 6) ' ) ' 

*DA3 * 

31 { 1 , : ) 

= 


(Bl(l, I), 1=1, M) 

WRITE (*, 

' A3C , SE13.6) ' ) ' 

* CA3 * 

31(2, : ) 

= 


(B1 (2,1) , 1 = 1, M) 

WRITE (*, 

' : A3 0 , 6E1 3 .6) ' ) ' 

* CA3 * 

31(3, : ) 

= 


(B1 (3,1), 1=1, M) 

WRITE (*, 

• >A30 , 613) 1 ) 

* DA 3 * 

UMAX 

= 


(UMAX ( I ) , 1=1, M) 

WRITE (*, 

' ■ A30, 6F14 ' 

*DA3 * 

XUMAX 

= 


( XUMAX ( I ) , 1=1, M) 

WRITE (*, 

' ( A3 0 , 318 . 6 > ’ ) 

*DA3 * 

TOLANG 

- 


TOLANG 

WRITE (*, 

' : A3 C , 3 3 1 3 . 6 ! ’ ) 1 

*DA3 * 

MXMAX 



( MXMAX ( I ) , 1=1,3) 

WRITE (*, 

' ■ A3 C , E13 . 6) ' ) 

* DA3 * 

MAXNORM 

= 


MAXNORM 

WRITE (*, 

' ! A 3 0 , E 1 3 . 6 i ’ ) 1 * 

DA3 * TC 

>L 

/ 

TOL 

WRITE (*, 

' i'/A50 ) 1 ) ' First: 

call t 

o R z 0 1 





END IF 

CALL R20 (Ul, IU, INFRONT , TEMP2 , TEMP3 , TEMP4 , ISVERTEX, 
Sc Bl, UMAX, XUMAX, U_MIN,U_MAX, TOL, M, DIAGS) 


(DIAGS) 

THEN 




WRITE (*, 

. ’ i/ASC)') ' 

A f 2. ci r i 3u K-i 1 



WRITE (*, 

, ’ ! A3 0 , 7I3i ' ] 

1 ' * DA3 * TEMP 2 ( I THETA) = 

' , (TEMP2 (I) , 

, 1=1 , TEMP3+1 

WRITE ( * , 

, ’ ;A3Q, 13) 1 ) 

1 *DA3* TEMP 3 i KANGS ) 

' , TEMP3 


WRITE ( * , 

, ’ ; A3 0 , 13) 1 ) 

( *DA3 * TEMP 4 ilMDX) 

’ , TEMP4 



END IF 

IF (IU.NE.0) THEN 

IF ( INFRONT . EQ . 1 . ) THEN 
DO 1=1, M 

U_L I F ( I ) = U1(I) 

ENDDO 

I_LIF = IU 

ELSEIF ( INFRONT. EQ. -1 . ) THEN 
DO 1=1, M 

U_L I B ( I ) = U1(I) 

ENDDO 

I_LIB = IU 
END IF 
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0286 


IF (DIAGS) THEN 


0287 


WRITE (*, 

' ;/A5C) ' ) ' After 1st L:F/L23' 


0288 


WRITE (*, 

•; A3 0,613)') ' *DA $* U L..F = (U_ 

LIF(I), 1= 

0289 


WRITE (*, 

■ ;A3C, i3) • ) ■ *;:a3* l:f - i_lif 

0290 


WRITE (*, 

• fA3C,6I3) ‘ ) 1 *CA3 * UJ...I3 - (U_ 

LIB(I), 1= 

0291 


WRITE (*, 

1 • A 3 C , 1 3 ! 1 ) ' *CA >* I LI3 * I LIB 

0292 


END IF 



0293 


END IF 



0294 

C 




0295 


IF ( ISVERTEX) THEN 


0296 

c 

WRITE (*, 

*) 1 TIME = TIME, 1 FIRST CALL i 

0 R2 0 ' 

0297 


CALL DOVERT(UDA, SAT, Ul, B , U_MIN, U_MAX, M, NORM) 


0298 


RETURN 



0299 


END IF 



0300 

c 




0301 

c 




0302 

c 

1st rotation about x-axis 


0303 

c 




0304 


IF (M.GE.8) 

THEN 


0305 


I COUNT = 

1 


0306 


ELSE 



0307 


I COUNT = 

2 


0308 


END IF 



0309 


I COUNT = 

1 


0310 


COSPHI = CSPHI (1, ICOUNT) 


0311 


SINPHI = INFRONT*CSPHI (2, ICOUNT) 


0312 


T22 ( 1 , 1 ) = 

COSPHI 


0313 


T22 ( 1 , 2 ) = 

-SINPHI 


0314 


T22 ( 2 , 1 ) = 

SINPHI 


0315 


T22 ( 2 , 2 ) = 

COSPHI 


0316 


MAXSTEPS = 

2 * INT ( ABS ( PI /ASIN (SINPHI) ) ) 


0317 


WAS INF = INFRONT 


0318 


ISOK = .FALSE. 


0319 


NMAX = NBI 

+ 1 


0320 


DIDSWITCH = 

. FALSE . 


0321 


STEPS = 0 



0322 


STUCK = .FALSE. 


0323 

c 




0324 


IF (DIAGS) 

THEN 


0325 


WRITE (*, 

' , ASG) ’ ) ' Before Main Loop' 


0326 


WRITE ( * , 

1 (A30.E1S.6) ') ' *SA3* COSPHI = ' 

, COSPHI 

0327 


WRITE ( * , 

' ( A3G, 318 . 6 S 1 ) • *DA3* SINPHI = 1 

, SINPHI 

0328 


WRITE (*, 

' ( A3 C , 1 3 ) ' ) 1 *DA3 * MAXSTEPS = ' 

, MAXSTEPS 

0329 


END IF 



0330 

c 




0331 

c 

MAIN LOOP ********★***********+*************** 


0332 

c 




0333 


DO WHILE ( ( ICOUNT . LT .NMAX) . AND . ( . NOT . ISOK) ) 


0334 

c 




0335 


STEPS = STEPS+1 


0336 


IF (STEPS 

.GE. MAXSTEPS) THEN 


0337 


STUCK = 

. TRUE . 


0338 

c 

WRITE {*, 

1 A3 0 1 1 ) * ★■*■★★*■****#**■*■*■# 1 


0339 

c 

WRITE (*, 

' i A3 0 , L3 ; ' ) ’ *CA3 * STUCK 

1 , STUCK 

0340 

c 

WRITE (* , 

' (A30, F14.6; ' ) ' *DA3* TIME 

' , TIME 

0341 

c 

WRITE ( * , 

' T T |'^ / , » J > *★**■**■*»*★★★*★* 1 


0342 


END IF 
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0343 


0344 

DO J = 1 , M 



0345 

BTEMP (1) 

= T22 (1,1) *B1 (2 , J) + T22 ( 1 , 2 ) 

*B1 (3, J) 

0346 

BTEMP (2) 

= T22 (2, 1) *B1 (2, J) + T22 (2,2) 

*B1 (3, J) 

0347 

B1 (2, J) 

= BTEMP (1) 


0348 

Bl(3, J) 

= BTEMP (2) 


0349 

ENDDO 



0350 C 




0351 

IF (DIAGS) 

THEN 


0352 

WRITE (*, 1 

t/A50) ' ) • In CA3 DOWHILE 1 


0353 

WRITE (*, ’ 

A3 0, 1 3 11 ) ' »DA3* I COUNT 

= ' , ICOUNT 

0354 

WRITE (*, ’ 

A3 C ,13! ' ) ' *DA3* STEPS 

= 1 , STEPS 

0355 

WRITE ( * , ' 

■A3C,I3> ') ' *DA3* XAXSTE? 

S = ' , MAXSTEPS 

0356 

WRITE (*, 1 

A3 0 , 6 E 1 8 . 6 ? 1 ) ' * DA3 * 31(2,:). 

= (B1 ( 2 , I ) 

0357 

WRITE (*, ' : 

A 1 0 i 6 — _ 3 . ) 1 ) ' * S A j B 1 ■ 3 i • ) 

= (B1 (3,1) 

0358 

ENDIF 




0359 

0360 

0361 

0362 

0363 

0364 

0365 

0366 

0367 

0368 

0369 

0370 

0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 

0386 

0387 

0388 

0389 

0390 

0391 

0392 

0393 

0394 

0395 

0396 

0397 

0398 

0399 


NANGS = TEMP3 
INDX = TEMP4 
DO 1=1,21 

ITHETA ( I } = TEMP2 (I) 

ENDDO 

CALL R20 (Ul, IU, INFRONT, TEMP2 , TEMP3 , TEMP4 , ISVERTEX, 
8c B 1 , UMAX , XUMAX , U MIN,U MAX, TOL, M, DIAGS) 


IF (DIAGS) THEN 
WRITE (*, ' : /A5'C) * ) 
WRITE <*, ' :A30, 713 
WRITE (*, r ; A3 0 , 13) 
WRITE (*, ' : A30 , 13) 
END IF 


\f.zer Loc 
' *DA3 * 

' *DA3* 

' *DA3* 


R20 1 
TEMPS 
TEMP 3 
TEMP 4 


;NAN G6 
( TNDX ) 


(TEMP2 (I) , 1= 

TEMP3 

TEMP4 


IF (IU.NE.0) THEN 

I F ( INFRONT . EQ . 1 . ) THEN 
DO 1=1 , M 

U_L I F { I ) = Ul (I) 

ENDDO 

I_LIF = IU 

ELSEIF ( INFRONT. EQ. -1. ) THEN 
DO 1=1, M 

U_LIB(I) = U1(I) 

ENDDO 

I_LIB = IU 
ENDIF 


CF (DIAGS) 
WRITE (*, 1 

THEN 

' After Loop R20 LIE/ 

L 13 1 


WRITE (*, 1 

i A3 0 , 6 1 3 • 

’) ' *CA3 * 'J _ DIF - 

(U_LIF(I) , 

1=1, M) 

WRITE (* , ’ 

i A3 G , 13 i ' 

) ' *DA3* I L I F = ' , 

I_LIF 


WRITE (*, ' 

A3 C ,613 .! 

1 ) 1 »CA3* *J_LI3 « 1 , 

(U_LIB (I) , 

1=1, M) 

WRITE (*, ' 

A3 C , 13) ’ 

) 1 *DA3* I LIB = 1 , 

I LIB 



END IF 
ENDIF 

IF (ISVERTEX) THEN 

WRITE (*,*) ’ f I ME = ’ , TIME, ’ LOOP CALL TO : 

CALL DOVERT ( UDA , SAT , Ul , B 1 , U MIN,U MAX , M, NORM) 


1,TEMP3+1) 
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0400 

RETURN 




0401 

END IF 




0402 





0403 

IF (DIAGS) THEN 



0404 

WRITE (*, ' 

(/ABO) •) ’ Refer 

e cesr. i n g reve r s a 1 ' 


0405 

WRITE (*, ' 

( A3 0 , FL4 .6) ' ) 

* DA 3 * IN FRONT - 

INFRONT 

0406 

WRITE ( * , 1 

■:a3g, ? 1 4 .6* ' ) 

*DA3* WAS INF - ’ , 

WAS INF 

0407 

END IF 




0408 





0409 

DIDSWITCH = 

.FALSE. 



0410 





0411 

IF (INFRONT. 

NE. WAS INF) THEN 

1 REVERSE DIRECTION 


0412 

IF (DIAGS) 

THEN 



0413 

WRITE (* 

, • C./A5C) ’) ' Rev 

arsing 1 


0414 

WRITE ( * 

, ' A3C, 13 ! ' ) 

’ *CA3* Steps taken 

- ' , STEPS 

0415 

WRITE (* 

, ' ’’ A3 C , FI 6 . 8 ’ ) 

' * DA 3 * Angle 

= 180 . *ASIN(SINPHI) /PI 

0416 

WRITE ( * 

/ ‘ A3 G , - j ; 1 ) 

' *DA3* MAXSTEPS 

= ' , MAXSTEPS 

0417 

WRITE (* 

, ' • A3 0 , L 3 • 1 ) 

' * DA3 * STUCK 

= ’ , STUCK 

0418 

END IF 




0419 

DIDSWITCH 

= . TRUE . 



0420 

WAS INF = INFRONT 



0421 

I COUNT = I COUNT +1 



0422 C 

Bisection and next 

transformation 




0423 

COSPHI = CSPHI (1, ICOUNT) 



0424 

SINPHI = INFRONT*CS PHI (2, ICOUNT) 


0425 

T2 2 (1,1) = 

COSPHI 



0426 

T22 (1,2) = - 

SINPHI 



0427 

T22 (2,1) = 

SINPHI 



0428 

T22 (2,2) = 

COSPHI 



0429 

MAXSTEPS = 2 

* INT (ABS (PI/ASIN (SINPHI) ) ) 


0430 

IF (DIAGS) THEN 



0431 

WRITE { * , ' 

! /A50 ) ' ) ' Bis 

action and next t 

r a ns format i. o n. ’ 

0432 

WRITE ( * , • 

:A3C, El 3 . 6; 1 ) 

1 *DA3* COSPHI 

- 1 , COSPHI 

0433 

WRITE (*, ' 

A30, E13 .6 i ' ) 

1 + CA3 * SINPHI 

- ’ , SINPHI 

0434 

WRITE ( * , 1 

' A3C , IS) ' ) 

' *DA3* MAXSTEPS 

- ’ , MAXSTEPS 

0435 

END IF 




0436 

STEPS = 0 




0437 C Check 

last edge with 

new B1 



0438 

Y = 0 . 0 




0439 

DO I = 1,M 




0440 

Y = Y + B1 (2, I) *XUMAX(I) 



0441 

ENDDO 




0442 

SY = 1 




0443 

IF (Y.LT. 0.0) SY = -1 



0444 





0445 

IF (DIAGS) THEN 



0446 

WRITE (*, ' : 

Ab 0 j 1 ) ' Be to r 

e GETEDGE’ 


0447 

WRITE ( * , ' :a 

30, 6 El 8 .6; ' ) ' 

*DA3 * 31 f' 2 , : / - 

’ , (B1(2,I) , 1= 1 , M) 

0448 

WRITE (*, ' ; A 

3 0 , 6 E 1 3 . 6 ii ' ) f 

* DA 3 * 31(3,:) = 

’ , ( B1 ( 3 , I ) ,1=1, M) 

0449 

WRITE (*, 1 (A 

3 0 , 6 l 3 .! r ) ’ 

*CA3 * UMAX 

’ , (UMAX ( I ) , 1=1, M) 

0450 

WRITE (*, • (A 

30, £13.6 - 1 ) 1 

^ ^ * ■/ ... 

' , Y 

0451 

WRITE!*, ' A 

30,13) 1 ) 

* D e 3 * 3 v ■' 

’ , SY 

0452 

WRITE (*, ■ ,:a 

3 0 , "i 1 3 i ' ) 

*DA3 * I THETA * 

' , ( ITHETA ( I ) , 1 = 1 , NANGS+1 

0453 

WRITE (*, ' , A 

30,13!’) 

*DA3* NANG 3 

* , NANGS 

0454 

END IF 





0455 

0456 CALL GETEDGE (U2 , JU, INFRONT, ISVERTEX, 
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0457 

0458 

0459 

0460 

0461 

0462 

0463 

0464 

0465 

0466 

0467 

0468 

0469 

0470 

0471 

0472 

0473 

0474 

0475 

0476 

0477 

0478 

0479 

0480 

0481 

0482 

0483 

0484 

0485 

0486 

0487 

0488 

0489 

0490 

0491 

0492 

0493 

0494 

0495 

0496 

0497 

0498 

0499 

0500 

0501 

0502 

0503 

0504 

0505 

0506 

0507 

0508 

0509 

0510 

0511 

0512 

0513 


C 


c 

c 


Sc Bl, UMAX, Y, SY, ITHETA, NANGS , U_MIN, U_MAX, INDX, TOL, M, DIAGS) 


IF (JU.NE.O) THEN 

IF ( INFRONT. EQ. 1. ) THEN 
DO 1=1, M 

U_L I F ( I ) = U2 (I) 

ENDDO 

I_LIF = JU 

ELSEIF ( INFRONT. EQ. -1 . ) THEN 
DO 1=1, M 

U_LIB (I) = U2 (I) 

ENDDO 

I_LIB = JU 
END IF 

IF (DIAGS) THEN 

WRITE (*, ‘ /A5G) ' ) ' After GETEDGE 

WRITE (*, ‘ A3C, 613) ' ) 1 *DA3* U_LIF - 

WRITE (*, 1 :A3C, 13) ' ) 1 + DA3 + I_LIF = 

WRITE ( * , ' . A3 C , 6 I 3 ) ') ' *DA3* U_LI3 = 
WRITE ( * , ' : A3 0 , 1 3 ) ' ) 1 *DA3* I_LI3 = 
END IF 
END IF 


( U_L I F ( I ) , 1=1, M) 
I_LIF 

(U_LIB(I), 1=1, M) 
I LIB 


IF ( ISVERTEX) THEN 

WRITE (* , * ) ' TIME = ' , TIME, ’ FOL GETEDGE 1 

CALL DO VE RT ( UD A , S AT , U2 , B , U_M I N , U_MAX , M , NORM ) 
RETURN 
END IF 


IF (JU.NE.O) THEN 

CALL IS FACET ( ISOK, IUOUT, JUOUT, U123, 

& IU, JU, Ul, U2 , Bl, M, TOL) 

IF (.NOT. ISOK) CALL ISFACET ( ISOK, IUOUT, JUOUT, U123, 
Sc I_LIF, I_LIB, U_L I F , U_LIB, Bl, M, TOL) 

ELSE 

ISOK = .FALSE. 

END IF 


END IF i IF (INFRONT. NS. WAS INF) THEN 


C 

C Must leave on a switch 
C 

C IF ( (ICOUNT.EQ.NMAX) .AND. ( .NOT. ISOK) .AND. ( .NOT. DIDSWITCH) ) THEN 

C NMAX = NMAX+ 1 

C END IF 


C 

C 


C 


IF (STUCK) THEN 

WRITE ( * , ' ■; / A5 0 } ’ ) 1 Stuck in DA 3 , exiting' 

WRITE (*,’•; A3 0 , FI 4 . 6 f ' ) ■ T I ME = ' , TIME 

RETURN 
END IF 

ENDDO ! End of do while statement 


C END MAIN LOOP *********+********************** 

C 
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0514 



IF (DIAGS) THEN 

0515 



WRITE (*, ' :/A50! ‘ ) ’ Zxicea frcr CA3 ' 

0516 



END IF 

0517 




0518 



IF ( . NOT . ISOK) CALL ISFACET ( ISOK, IUOUT, JUOUT, U123, 

0519 


Sc 

I LIF, I LIB , U LIF, U LIB, Bl, M, TOL) 

0520 



IF (ISOK) THEN 

0521 



DO 1=1,3 

0522 



DO J=1 , M 

0523 



IUTEMP (J) =U12 3 (J, I) 

0524 



ENDDO 

0525 



CALL SETU(XUTEMP, IUTEMP, U MIN,U MAX, M) 

0526 



DO J=1 , M 

0527 



XU123 (J, I) =XUTEMP ( J) 

0528 



ENDDO 

0529 



ENDDO 

0530 

C 



0531 



DO 1=1,3 

0532 



DO J=1 , 3 

0533 



M123 (I, J) =0 . 

0534 



DO K=1 , M 

0535 



M123 (I, J)=M123 (I, J) +B(I,K) *XU123 (K, J) 

0536 



ENDDO 

0537 



ENDDO 

0538 



ENDDO 

0539 

c 



0540 



DO 1=1,3 

0541 



DO J=2 , 3 

0542 



M123 (I, J) =M123 (1,1) -M123 (I, J) 

0543 



ENDDO 

0544 



MTEMP (I) =M123 (1,1) 

0545 



M12 3 (1,1) =MDES (I) 

0546 



ENDDO 

0547 

c 



0548 



CALL INVMAT3 (M123 , MINV, DET) 

0549 

c 



0550 



DO 1=1,3 

0551 



ABC (I) = 0. 

0552 



DO J=l, 3 

0553 



ABC ( I ) = ABC ( I) +MINV ( I , J) *MTEMP ( J) 

0554 



ENDDO 

0555 



ENDDO 

0556 



AA = ABC ( 1 ) 

0557 



BB = ABC (2) 

0558 



CC = ABC (3) 

0559 



SAT = l./AA 

0560 



IF (AA.LT. 1 . ) AA = 1 . 

0561 



DO 1=1, M 

0562 



UDA(I) = (XU123 (1,1) 

0563 


Sc 

+BB* (XU123 (1,2) -XU123 (1,1)) 

0564 


Sc 

+CC* (XU123 (1,3) -XU123 (1,1))) /AA 

0565 



ENDDO 

0566 




0567 



IERR = 0 

0568 

c 



0569 

c 

Call 

estimate subroutine to estimate solution if facet not 

0570 

c 




found 
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0571 

0572 

0573 

0574 

0575 

0576 

0577 

0578 

0579 

0580 

0581 

0582 

0583 

0584 

0585 

0586 

0587 

0588 

0589 

0590 

0591 

0592 

0593 

0594 

0595 

0596 

0597 

0598 

0599 

0600 
0601 
0602 

0603 

0604 

0605 

0606 

0607 

0608 

0609 

0610 
0611 
0612 

0613 

0614 

0615 

0616 

0617 

0618 

0619 

0620 
0621 
0622 

0623 

0624 

0625 

0626 
0627 


ELSE 

CALL EST (UDA, SAT, I ERR, 

& U_L I F , I_LIF, U_L I B , I_LIB, Bl, U_MIN, U_MAX, NORM, M) 
END IF 


RETURN 

END 


SUBROUTINE EST (UDA, SAT, IERR , 

Sc Ul, IU, U2 , JU, B, U_MIN, U_MAX, NORM, M) 
IMPLICIT NONE 


INPUTS 

REAL *4 B ( 3 , 2 0 ) , U_MAX(20), U_MIN(20), NORM 
INTEGERM Ul(20), IU, U2(20), JU, M 

OUTPUTS 

REALM SAT, UDA(20) 

INTEGER* 4 IERR 


LOCALS 

REALM XU1(20), XU2 (20) , XU3 (20) , XU4 (20) , XMOM(3) 

REAL*4 UPPER1 ( 3 ) , UPPER2(3), LOWER1 ( 3 ) , LOWER2(3), XNORM 
REAL*4 XK1 , XK2 , XK3 , XVI (3), XV2(3), XW1 (20) , XW2(20) 
INTEGER*4 U3 (20) , U4(20) 


OTHER LOCALS 

INTEGER*4 I, J, K 


IERR = 1 


Ul(IU) = -1 
U2 ( JU) = -1 
DO 1=1 , M 

U3 ( I ) = Ul ( I ) 
U4 (I) = U2 (I) 
ENDDO 

U3(IU) = 1 

U4 ( JU) = 1 


CALL SETU ( XU1 , Ul , U_MIN , U___MAX , M) 
CALL SETU (XU2 , U2 , U_MIN, U_MAX , M ) 
CALL SETU (XU3 , U3 , U_MIN, U_MAX , M) 
CALL SETU (XU4 , U4 , U MIN,U MAX , M) 


DO 1=1,3 

LOWER1 ( I ) = 

LOWER2 ( I ) = 
UPPERl(I) = 
UPPER2 ( I ) = 

DO J= 1 , M 
LOWER1 (I) 
LOWER2 (I) 
UPPER1 (I) 
UPPER2 (I) 


0 . 
0 . 
0 . 
0 . 


= LOWER1 (I) +B(I, J) *XU1 (J) 
= LOWER2 (I) +B(I, J) *XU2 (J) 
= UPPER1 (I) +B (I, J) *XU3 (J) 
= UPPER2 (I) +B ( I , J) *XU4 ( J) 



Alloc . f 

Printed: Friday, March 23, 2001 12:46:16 PM 


Page 12 of 28 
Printed For: Bull Durham 


0628 

0629 

0630 

0631 

0632 

0633 

0634 

0635 

0636 

0637 

0638 

0639 

0640 

0641 

0642 

0643 

0644 

0645 

0646 

0647 

0648 

0649 

0650 

0651 

0652 

0653 

0654 

0655 

0656 

0657 

0658 

0659 

0660 
0661 
0662 

0663 

0664 

0665 

0666 

0667 

0668 

0669 

0670 

0671 

0672 

0673 

0674 

0675 

0676 

0677 

0678 

0679 

0680 
0681 
0682 

0683 

0684 


ENDDO 

ENDDO 

IF ( LOWER1 ( 2 ) . NE . UPPER1 ( 2 ) ) THEN 

XK1 = LOWER 1 (2 ) / (LOWER1 (2 ) -UP PERI ( 2 ) ) 
ELSE 

XK1 = 0 . 

END IF 

IF (LOWER2 (2) .NE.UPPER2 (2) ) THEN 

XK2 = LOWER2 (2) / (LOWER2 (2) -UPPER2 (2) ) 
ELSE 

XK2 = 0. 

ENDIF 


DO 1=1,3 

XVI (I) = XK1*UPPER1 (I) + (1 . -XK1) *L0WER1 (I) 

XV2 ( I ) = XK2*UPPER2 (I) + (1 . -XK2) *LOWER2 (I) 
ENDDO 

IF (XV2 (3) .NE.XV1 (3) ) THEN 
XK3 = XV2 (3) / {XV2 (3) -XVI (3) ) 

ELSE 

XK3 = 0 . 

ENDIF 

DO 1=1, M 

XW1(I) = XK1*XU3 (I) + (1. -XK1) *XU1 (I) 

XW2 { I ) = XK2*XU4 (I) + (1. -XK2) *XU2 (I) 

ENDDO 

DO 1=1, M 

UDA(I) = XK3*XW1 (I) + (1 . -XK3) *XW2 (I) 

ENDDO 


DO 1=1, 3 
XMOM ( I ) = 0 . 

DO J=1 , M 

XMOM (I) = XMOM ( I ) + B ( I , J ) * UDA ( J ) 

ENDDO 
ENDDO 

XNORM = SQRT (XMOM ( 1 ) *XMOM ( 1 } 

Sc +XMOM { 2 ) *XMOM ( 2 ) 

Sc +XMOM ( 3 ) *XMOM ( 3 ) ) 

IF (XNORM. NE. 0. ) THEN 
SAT = NORM/XNORM 
XNORM = SAT 
ELSE 

SAT = 0. 

ENDIF 

IF (XNORM. GT. 1 . ) XNORM = 1. 

DO I = 1 , M 
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0685 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 

0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 

0711 

0712 

0713 

0714 

0715 

0716 

0717 

0718 

0719 

0720 

0721 

0722 

0723 

0724 

0725 

0726 

0727 

0728 

0729 

0730 

0731 

0732 

0733 

0734 

0735 

0736 

0737 

0738 

0739 

0740 

0741 


UDA ( I ) = XNORM* UDA ( I ) 

ENDDO 

RETURN 

END 


SUBROUTINE DOVERT (UDA, SAT, 

Sc Ul, B,U_MIN,U_MAX,M,NORM) 

IMPLICIT NONE 

REALM UDA (20), SAT, U_MIN(20), U_MAX(20), B(3,20), NORM 
INTEGERM Ul(20), M 

REALM XMOM ( 3 ) , XNORM 
INTEGERM I, J 

WRITE (*,*) ' VERTEX ’ 

CALL SETU (UDA , U1 , U_MIN , U_MAX , M ) 

DO 1=1,3 

XMOM ( I ) = 0 . 

DO J= 1 , M 

XMOM (I) = XMOM (I) +B(I, J) *UDA(J) 

ENDDO 

ENDDO 


XNORM = SQRT (XMOM (1) *XMOM ( 1) 

Sc +XMOM ( 2 ) *XMOM ( 2 ) 

Sc +XMOM ( 3 ) *XMOM ( 3 ) ) 

SAT = NORM/ XNORM 
XNORM = SAT 

IF (XNORM. GT. 1. ) XNORM = 1. 

DO I = 1 , M 

UDA (I) * XNORM*UDA ( I ) 

ENDDO 

RETURN 

END 


SUBROUTINE INVMAT3 (MATIN, MATOUT, DET) 

IMPLICIT NONE 
INTEGERM I, J 

REALM DET, MATIN(3,3), MATOUT ( 3 , 3 ) 


Zero out the output matrix 

DO I = 1,3 
DO J = 1,3 
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0742 

0743 

0744 

0745 

0746 

0747 

0748 

0749 

0750 

0751 

0752 

0753 

0754 

0755 

0756 

0757 

0758 

0759 

0760 

0761 

0762 

0763 

0764 

0765 

0766 

0767 

0768 

0769 

0770 

0771 

0772 

0773 

0774 

0775 

0776 

0777 

0778 

0779 

0780 

0781 

0782 

0783 

0784 

0785 

0786 

0787 

0788 

0789 

0790 

0791 

0792 

0793 

0794 

0795 

0796 

0797 

0798 


MATOUT ( I , J) =0.0 
ENDDO 
ENDDO 
C 

C Calculate the determinant of the input matrix 
C 

DET = MATIN (1,1) * MAT IN (2,2) * MAT IN (3,3) 

5c + MATIN (1,2) * MATIN (2,3) * MATIN (3,1) 

& + MATIN (1,3) * MAT IN (2,1) * MATIN (3,2) 

5c - MATIN (1,3) * MATIN (2,2) * MATIN (3,1) 

& - MATIN (1,2) * MATIN (2,1) * MATIN (3,3) 

Sc - MATIN (1,1) * MATIN (2,3) * MATIN (3,2) 

C 

C Find the matrix inverse 


C 


IF (DET. NE. 0.0) 
MATOUT (1,1) = 
MATOUT (1,2) = 
MATOUT (1,3) = 
MATOUT (2,1) = 
MATOUT (2,2) = 
MATOUT (2, 3) = 
MATOUT (3,1) = 
MATOUT (3,2) = 
MATOUT (3,3) = 
END IF 


THEN 

(MATIN (2, 2) * MATIN ( 3 , 3) 
- (MATIN (1, 2) * MATIN ( 3 , 3) 
(MATIN (1, 2) * MATIN (2,3) 
- (MATIN (2, 1) *MATIN ( 3 , 3) 
(MATIN (1, 1) * MATIN ( 3 , 3) 
- (MATIN (1, 1) * MAT IN (2,3) 
(MATIN (2, 1) * MATIN ( 3 , 2) 
- (MATIN (1, 1) * MAT I N ( 3 , 2 ) 
(MATIN (1, 1) * MAT IN (2,2) 


-MATIN (2, 3) * MATIN (3,2) 
-MATIN (1, 3) ♦MATIN (3,2) 
-MATIN (1, 3) * MAT IN (2,2) 
-MATIN (2, 3) * MAT IN ( 3 , 1) 
-MATIN (1, 3) * MATIN ( 3 , 1) 
-MATIN (1, 3) * MATIN ( 2 , 1) 
-MATIN (2,2) * MAT IN ( 3 , 1) 
-MATIN (1,2) * MAT I N ( 3 , 1) 
-MATIN (1,2) * MAT IN ( 2 , 1) 


) /DET 
) /DET 
) /DET 
) /DET 
) /DET 
) /DET 
) /DET 
) /DET 
) /DET 


RETURN 

END 


C 

SUBROUTINE DCGEN (T, MD) 

IMPLICIT NONE 
REALM MD ( 3 ) 

REAL* 4 T ( 3 , 3 ) 

INTEGER* 4 MLOCAL 

REAL* 8 V ( 3 ) , VLEN ( 3 ) , VNORM, XDOT_DC 
REAL *4 DETNUM, AMIN_DC 

INTEGERM J, JCOL, KCOL, I, MOM_FLAG , IROW 
REAL* 8 DTOL 
C 

C Calculate the norm of the moments 
C 

DTOL = l.D-5 

VNORM = DSQRT (dble(MD(l) ) *dble (MD ( 1 ) ) 

& +dble (MD (2) ) *dble (MD (2) ) 

Sc +dble (MD ( 3 ) ) *dble(MD(3) ) ) 

C 

C Normalize the desired moments 
C 

DO I = 1,3 

V (I ) = dble (MD ( I ) ) /VNORM 
ENDDO 
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0799 

0800 
0801 
0802 

0803 

0804 

0805 

0806 

0807 

0808 

0809 

0810 
0811 
0812 

0813 

0814 

0815 

0816 

0817 

0818 

0819 

0820 
0821 
0822 

0823 

0824 

0825 

0826 

0827 

0828 

0829 

0830 

0831 

0832 

0833 

0834 

0835 

0836 

0837 

0838 

0839 

0840 

0841 

0842 

0843 

0844 

0845 

0846 

0847 

0848 

0849 

0850 

0851 

0852 

0853 

0854 

0855 


C 

C Zero out the transformation matrix 
C 

DO I = 1,3 
DO J = 1,3 
T ( I , J) =0.0 
ENDDO 
ENDDO 
C 

C Check to see if V(I),3 to 1 is approx equal to zero => reduce size of problem 
C 

DO I = 3, 1, -1 

IF (ABS (V { I ) ) .LE. DTOL) THEN 
T(I, I) = 1.0 
ELSE 
GOTO 5 
END IF 
ENDDO 
C 

C T ( 1 , 1 ) = 1.0 or -1.0 for rotation about x-axis (depends on direction of rot.) 
C 

5 IF (I .EQ. 1) THEN 
T ( I , I ) = 1.0 

IF (dble(MD(l) ) .LT.0.0D0) T(I,I) = -1.0 
RETURN 
END IF 
C 

C Set the 1st row of T equal to the normalized desired moments and 
C calculate the square of each of these values 
C 

MLOCAL = I 
DO I = 1,3 

T ( 1 , I ) = V (I) 

VLEN(I) = V ( I ) *V(I) 

ENDDO 

C 

C Developing orthogonal tranf ormation with V as 1st row 
C 

DO JCOL = 1 , MLOCAL- 1 

I ROW = MLOCAL + 1 - JCOL 

T ( I ROW, JCOL) = sngl (DSQRT ( 1 . 0D0 - VLEN(JCOL))) 

DO KCOL = JCOL+1, MLOCAL 

XDOT_DC = T(l, JCOL) *T(1, KCOL) 

IF (I ROW .NE. MLOCAL) THEN 
DO I = IROW+1 , MLOCAL 

XDOT_DC = XDOT_DC + T { I , JCOL) *T ( I , KCOL ) 

ENDDO 

ENDIF 

T { IROW, KCOL) = -XDOT_DC/T ( I ROW, JCOL) 

VLEN(KCOL) = VLEN(KCOL) + dble (T ( IROW, KCOL) ) *dble (T ( IROW , KCOL) ) 

ENDDO 

ENDDO 

C 

C Tricky stuff here! 

C 

DETNUM = int (mod (MLOCAL, 4) /2.Q) 
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0856 

0857 

0858 

0859 

0860 
0861 
0862 

0863 

0864 

0865 

0866 

0867 

0868 

0869 

0870 

0871 

0872 

0873 

0874 

0875 

0876 

0877 

0878 

0879 

0880 
0881 
0882 

0883 

0884 

0885 

0886 

0887 

0888 

0889 

0890 

0891 

0892 

0893 

0894 

0895 

0896 C 

0897 

0898 

0899 

0900 

0901 

0902 

0903 

0904 

0905 

0906 

0907 

0908 

0909 

0910 

0911 

0912 C 


" Mecessary to do, but not easy to explain 

IF (DETNUM . EQ . 0 . 0 ) 

IF ( T ( 1 , MLOCAL) 

T ( 2 , MLOCAL - 1 ) 

T (2, MLOCAL) 

END IF 
ELSE 

IF (T(l, MLOCAL) 

T (2, MLOCAL- 1) 

T( 2, MLOCAL) 

END IF 
END IF 


THEN 

.LT. 0.0) THEN 
= -T (2 , MLOCAL- 1 ) 
= -T (2, MLOCAL) 


.GT. 0.0) THEN 
= -T (2, MLOCAL- 1) 
= -T( 2, MLOCAL) 


RETURN 

END 


SUBROUTINE SETU (XU_SETU, IU_SETU, U_MIN, U_MAX , M) 

IMPLICIT NONE 

REAL *4 U_MAX (20), U_MIN(20) 

INTEGERM IMODE, M 
INTEGER* 4 I, IU_SETU(20) 

REAL* 4 XU_SETU (20) 

DO I = 1 , M 

IF ( IU_SETU ( I ) .EQ. 1) THEN 
XU_SETU(I) = U_MAX ( I ) 

ELSEIF ( IU_SETU { I ) . EQ . - 1 ) THEN 

XU_SETU(I) = U_M I N { I ) 

ELSE 

XU_SETU ( I ) = 0 . 

END IF 
ENDDO 

RETURN 

END 


SUBROUTINE R20 {Ul, IU, INFRONT, ITHETA, NANGS , INDX, ISVERTEX, 
Sc B1 , UMAX , XUMAX , U_MIN, U_MAX , TOL , M , DIAGS ) 

IMPLICIT NONE 

REAL*4 U_MIN (20), U_MAX(20) 

INTEGER*4 IMODE, ITHETA (21), IU, UMAX(20), M, N, NANGS , SY 
INTEGER*4 Ul(20) 

REAL*4 PIOVR2 , Bl(3,20), PI, XU(20), Y, XUMAX(20) 

REAL *4 INFRONT, TOL 
REAL *4 THETA (21) 

LOG I CAL* 4 ISVERTEX, DIAGS 
INTEGER*4 INDX, SSY 
INTEGER*4 I, J 

REAL* 4 ANG, YPREV, XU1(20), K, Z, Zl, Z2 
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0913 

0914 

0915 

0916 

0917 

0918 

0919 

0920 

0921 

0922 

0923 

0924 

0925 

0926 

0927 

0928 

0929 

0930 

0931 

0932 

0933 

0934 

0935 

0936 

0937 

0938 

0939 

0940 

0941 

0942 

0943 

0944 

0945 

0946 

0947 

0948 

0949 

0950 

0951 

0952 

0953 

0954 

0955 

0956 

0957 

0958 

0959 

0960 

0961 

0962 

0963 

0964 

0965 

0966 

0967 

0968 

0969 


DATA PIOVR2/1 . 570796326794897/ 


DATA PI/3 . 

141592653589793/ 




IF (DIAGS) 

THEN 




WRITE (* 

, ’ ; /A5C ) 1 ) 1 Enter 

ing R2C 



WRITE (* 

, ’ -A5C)') 4 Cal line 

g arguments' 



WRITE (* 

, ’ A 3 0 , 6E1 3 . 6 ) ‘ ) 4 

* p 2 o * 3i ■ ' ; 

= (Bid, i), 

1 = 1 , M) 

WRITE (* 

, 1 : A30, 6S1S . 6 } ' ) ' 

★ R 2 0 * 3 i ■ 2 : ) 

= (b: ( 2 ,d, 

1 = 1 , M) 

WRITE (* 

, 1 A 3 0 , 6 Z 1 3 .Si') 1 

* p 2 C * 3 1 ; 5 : i 

- ', (B1 (3 , I ) , 

1 = 1 , M) 

WRITE (* 

, 1 (A3C, 613! ' ) 

*R2C* Lie AX 

- (UMAX (I), 

1 = 1 , M) 

WRITE (* 

, ' : A3C , 6 ? It . 6 ) 1 ) ' 

*R2Q* XUMAX 

= ' , (XUMAX (I) 

, 1=1, M 

WRITE (* 

, ’ .'A3C , E18 . 6 i 1 ) 1 

★ P' ^ C * 

= ' , TOL 


WRITE (* 

, ' ■ AiO f -.3 1 ) 1 

★ p 2 C ■* v 

- ' , M 



ENDIF 


Calculate Y 

Y = 0 . 0 
DO I = 1 , M 

U1 (I) = UMAX (I) 

Y = Y + B1 (2 , I ) *XUMAX ( I ) 

ENDDO 

VERTEX CHECK 

ISVERTEX = .FALSE. 

IF (ABS ( Y) . LT . TOL) THEN 

Y = 0. 

Z = 0. 

DO 1=1, M 

Z = Z + Bl (3,1) * XUMAX ( I ) 

ENDDO 

IF (ABS (Z) .LT.TOL) THEN 

WRITE (*,*} ' Z ::: ’ ,Z, ’ TOL - ' , TOL , ' P.20, ABS £ Z ; . LT . TOL 

ISVERTEX = .TRUE. 

IU = 0 

INFRONT = 0 . 

NANGS = 0 
RETURN 
ENDIF 
ENDIF 

END VERTEX CHECK 

IF (ABS (Y) .LT.TOL) THEN 
SY = 0 

ELSEIF ( Y . LT .0.0) THEN 
SY = -1 
ELSE 

SY = 1 

ENDIF 

IF (DIAGS) THEN 

WRITE { * , ' ( . A5 3 j ' ) ' First calculations ’ 

WRITE (*, ’ ;A30, £13 .6; ' ) ’ *R2G* Y = ' , Y 

WRITE ( * , ’ ; A3 0 , 3 ’ ) ' *R20* SY = 1 , SY 

ENDIF 


Get the angle 
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0970 

C 


0971 

c 


0972 

IF (DIAGS) THEN 


0973 

WRITE (*,' /A5C) ’ ) 1 Gearing angles' 


0974 

END IF 


0975 

C 


0976 

NANGS = 0 


0977 

DO 1=1, M 


0978 

ANG = ATAN2 (-B1 (1, I) ,B1 (2, I) ) 


0979 

IF (ABS (ANG) .GT. PIOVR2) THEN 


0980 

IF (ANG.LT.O.) THEN 


0981 

ANG = ANG+PI 


0982 

ELSE 


0983 

ANG = ANG- PI 


0984 

END IF 


0985 

ENDIF 


0986 

NANGS = NANGS -hI 


0987 

THETA (NANGS) = ANG 


0988 

I THETA (NANGS) = I 


0989 

IF (DIAGS) THEN 


0990 

WRITE (*, ' :A30, 13, A10, E13 . 6 ) ' ) 1 *R20* 1 - 

’ / I / ' AN 

0991 

ENDIF 


0992 

ENDDO 


0993 

C 


0994 

0995 

C THETA and ITHETA now sorted by control number 


0996 

C Sort THETA by magnitude. ITHETA gets shuffled the 

same way 

0997 

C 


0998 

IF (DIAGS) THEN 


0999 

WRITE (*, ’ ; /AS 0 ) ') ' Before sorting' 


1000 

WRITE (* , ' : A30 , 13 ) 1 ) ' *R20* NANGS = ! , 

NANGS 

1001 

WRITE (*,' :A30, 613) ' ) ’ *R2C* ITHETA - ', 

( ITHETA ( I ) , 

1002 

WRITE ( * , ' ( A3 0 , 6 E 1 3 . 6 ) ' ) ' *R20* THETA 

’ , (THETA (I) 

1003 

ENDIF 


1004 

CALL SORTC (THETA, ITHETA, NANGS, THETA, ITHETA) 


1005 

IF (DIAGS) THEN 


1006 

WRITE (*, ’ A50) 1 ) ' After sorting’ 


1007 

WRITE (*, 1 ( A3 0 , 13) ' ) ’ *R2G* NANGS = ' , 

NANGS 

1008 

WRITE (*,' (A30, 613) ' ) 1 *R.2C* ITHETA ■= 

(ITHETA (I) , 

1009 

WRITE (*, ’ iA30 f 6S13 . 6} ' ) ' *R20* THETA = 1 

(THETA (I) 

1010 

ENDIF 


1011 

C FIND INDEX OF ZERO 


1012 

C 


1013 

DO 1=1, NANGS 


1014 

IF (THETA (I) .GT. 0 . ) THEN 


1015 

INDX = I 


1016 

DO J=NANGS , 1,-1 


1017 

THETA (J+l) = THETA (J) 


1018 

ITHETA (J+l) = ITHETA (J) 


1019 

ENDDO 


1020 

THETA (I) = 0. 


1021 

ITHETA (I) = 0 


1022 

GOTO 193 


1023 

ENDIF 


1024 

ENDDO 


1025 

193 CONTINUE 


1026 

C 



- 1 , ANG 


=1 , NANGS) 
1=1, NANGS ) 


=1, NANGS) 
1=1, NANGS) 
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1027 IF (DIAGS) THEN 

1028 WRITE ( * , ' : / AS 0 ■ ' ) ' Afte- ; naex ing 1 

1029 WRITE ( * , 1 : A 3 C , I 3 ' ' ) ' *R2C* JMDX = 1 , INDX 

1030 WRITE ( * , * : A 3 C , 7 1 3 ! ' ) ' *F.2C* ITHETA - 1 , ( ITHETA ( I ) , I=1,NANGS+1) 

1031 WRITE (*, ' :A3C,7E1S.6) 1 ) ' *R2C* THETA - ', (THETA (I) , I=1,NANGS+1) 

1032 END IF 

1033 C 

1034 CALL GETEDGE (U1 , IU, INFRONT, ISVERTEX, 

1035 C AS FUNCTIONS OF 

1036 & Bl, UMAX, Y, SY, ITHETA, NANGS , U_MIN, U_MAX, INDX, TOL , M, DIAGS) 

1037 c 

1038 IF (DIAGS) THEN 

1039 WRITE (* , ' i /A50 ) ‘ ) 1 Exiling R2C 1 

1040 WRITE ( * , ‘ ■ A3 C , 6 1 2 • ' ) ' *R2C* U1 - ', (U1(I), 1 = 1, M) 

1041 WRITE (*, ' A3 0 , I 3 ■ ‘ ) 1 *R2C* _U = \ IU 

1042 WRITE ( * , 1 ■ A3 i , Eln . 6 1 ) 1 * ? 2 C * I N FRONT - 1 , INFRONT 

1043 WRITE (*, 1 ;A30, T3 ■• ) 1 *R2C* ISVEPTEX - ' , ISVERTEX 

1044 END IF 

1045 c 

1046 RETURN 

1047 END 

1048 

1049 C 

1050 

1051 SUBROUTINE GETEDGE (U2 , JU, INFRONT, ISVERTEX, 

1052 C AS FUNCTIONS OF 

1053 & Bl, UMAX, Y, SY, ITHETA, NANGS, U_MIN, U_MAX, INDX, TOL, M, DIAGS) 

1054 

1055 IMPLICIT NONE 

1056 

1057 INTEGERM M, SY, UMAX(20), INDX, NANGS, ITHETA ( 21 ) 

1058 REAL *4 Bl(3,20), U_MAX(20), U_MIN(20), XU(20), Y 

1059 REAL* 4 INFRONT, TOL 

1060 INTEGER*4 U2(20), JU, SSY 

1061 LOG I CAL *4 ISVERTEX, DIAGS, STUCK 

1062 INTEGER*4 I, J, ICOUNT, IX 

1063 REAL*4 XK, XU2(20), YPREV, Z, Zl, Z2 , SAVE_Y, MOM (3) 

1064 C 

1065 IF (DIAGS) THEN 

1066 WRITE (*,' i/ASC) 1 ) ' Entering GETEDGE • 

1067 END IF 

1068 C 

1069 DO 1 = 1, M 

1070 U2 (I) = UMAX ( I ) 

1071 ENDDO 

1072 C 

1073 SAVE_Y = Y 

1074 C 

1075 IX = INDX - SY 

1076 SSY = SY 

1077 C 

1078 STUCK = .FALSE. 

1079 ICOUNT = 0 

1080 C 

1081 IF (DIAGS) THEN 

1082 WRITE <*,’: AS A 1 ) ' Beginning GETEDGE DCWH IDE 1 

1083 WRITE ( * , r .. AT 3,613 ’ ) 1 * GST EDGE* GT - ’, (U2 ( I ) , 1 = 1, M) 



Alloc . f 

Printed: Friday, March 23, 2001 12:46:16 PM 


Page 20 of 28 
Printed For: Bull Durham 


1084 

1085 

1086 

1087 

1088 

1089 C 

1090 

1091 

1092 

1093 

1094 

1095 

1096 

1097 

1098 

1099 

1100 
1101 
1102 

1103 

1104 

1105 

1106 

1107 C 

1108 

1109 

1110 
1111 
1112 

1113 

1114 

1115 

1116 

1117 

1118 

1119 

1120 
1121 
1122 

1123 

1124 

1125 

1126 C 

1127 

1128 C 

1129 

1130 C 

1131 

1132 

1133 C 

1134 

1135 

1136 

1137 

1138 

1139 

1140 


WRITE ( * , ' : A 3 C , Z 1 3 . 6 ■ 1 ) 

WRITE (*, 1 • A 3 C , 13 ) ’ ) 

WRITE (*, 1 :A30, 13} ' ) 

WRITE (*, 1 ( A 3 0 , 13) ' ) 

END IF 

DO WHILE ( (SY.EQ.SSY) .AND. (IX.GT.0) .AND. { IX . LE . NANGS ) ) 
I COUNT = ICOUNT+ 1 
YPREV = Y 
JU = I THETA ( IX) 

Y = Y-U2 ( JU) *B1 (2 , JU) * (U_MAX { JU) -U_MIN ( JU) ) 

SSY =1 

IF (Y.LT.0.0) SSY = -1 
IF (ABS (Y) .LT.TOL) THEN 
Y = 0 . 

SSY = 0 
END IF 

U2 ( JU) = -U2 ( JU) 

IF (SY.NE.0.) THEN 
IX = IX-SY 
ELSE 

IX = IX-1 
END IF 


, SAVE_Y 
, IX 
, SSY 
, NANGS 


I COUNT = I COUNT +1 


IF (ICOUNT 

GT .20) STUCK = 

. TRUE . 



IF (DIAGS) 

THEN 




WRITE ( * , ’ 

/AS 05 ' ) ' Botto 

n of GETEDG 

E DC WHILE ' 


WRITE (*, 1 

A3C , 13) ' ) 

* GETEDGE* 

I COUNT 

- ' , ICOUNT 

WRITE (* , ’ 

A3C , L3) ' ) 

* GETEDGE* 

STUCK 

= ' , STUCK 

WRITE (*, ’ 

A30 , E18 . 6 ) ' ) 1 

* GETEDGE* 

YPREV 

= 1 , YPREV 

WRITE (*, ' 

A30 , 13) ' ) 

* GETEDGE* 

JU 

= 1 , JU 

WRITE (*, 1 ■ 

A30 , E13 . 6 ) ’ ) ' 

* GETEDGE * 

y 

- 1 , Y 

WRITE (*, ' 

A3 0 , 13} ' ) 

* GETEDGE* 

SSY 

= ’ , SSY 

WRITE (*, ' 

A3 C , L3 } ' ) 

* GETEDGE* 

S '{ . u . S S Y 

= 1 , (SY.EQ.SSY) 

WRITE (*, ' 

A3 C , 6 1 3 ! ' ) 

* GETEDGE* 

U2 

= 1 , (U2 (I) , 1=1, M) 

WRITE (*, ’ 

A3C , 13) ' ) 

* GETEDGE* 

’ V 

= 1 , ix 

END IF 
IF (STUCK) 

THEN 




PAUSE ( ' St 

uck in GETEDGE 

in COWHIDE ’ 

) 


STUCK = . 

FALSE . 





END IF 


ENDDO 

IF (SY.EQ.SSY) THEN 

STUCK = .FALSE. 

I COUNT = 0 


IF (DIAGS) THEN 
WRITE (*, ' :/A50) ’ ) 1 

WRITE (*, 1 ;A3u, 613 i ' ) 
WRITE ( * , 1 :A3C, E 1 8 . 6 i 
WRITE (*, ' • A 2 C , 13: 1 ) 
WRITE (*, ' : A3 0 , 13; 1 ) 
WRITE ( * , 1 : A3 3 , 13.' ’ ) 



*■ 3 £ T £ Oc‘- £ * ;.£ 


DOWH HE 1 

- \ (U2 ( I ) , 1=1, M) 

Y - r , SAVE_Y 

= ' , IX 

- 1 , SSY 

■= * , NANGS 
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1141 

1142 C 

1143 

1144 

1145 

1146 

1147 

1148 

1149 

1150 

1151 

1152 

1153 

1154 

1155 

1156 

1157 

1158 

1159 

1160 

1161 C 

1162 

1163 

1164 

1165 

1166 

1167 

1168 

1169 

1170 

1171 

1172 

1173 

1174 

1175 

1176 

1177 

1178 

1179 C 

1180 
1181 

1182 C 

1183 

1184 

1185 C 

1186 

1187 

1188 

1189 C 

1190 

1191 

1192 

1193 

1194 

1195 

1196 

1197 


ENDIF 

Y = SAVE_Y 
IX = INDX+SY 

DO WHILE ( {SY.EQ.SSY) .AND. (IX.GT.0) .AND. ( IX . LE . NANGS+ 1 ) ) 
YPREV = Y 
JU = I THETA ( IX) 

Y = Y-U2 ( JU) *B1 (2 , JU) * (U_MAX ( JU) -U_MIN ( JU) ) 

SSY = 1 

IF { Y . LT . 0 .0) SSY = -1 
IF (ABS (Y) .LT.TOL) THEN 
Y = 0. 

SSY * 0 
ENDIF 

U2 ( JU) = -U2 ( JU) 

IF (SY.NE.0.) THEN 
IX = IX+SY 
ELSE 

IX = IX+1 
ENDIF 

I COUNT = I COUNT +1 

IF { I COUNT . GT .20) STUCK = .TRUE. 

IF (DIAGS) THEN 


WRITE ( * , ' 

s / A5 C ) ') 

1 

Bet con 

of Alt. G 

EG' EDGE DOWHILE ’ 


WRITE (*, 1 

A3C , 13) 1 

) 

I 

*GETEDG2* 

~ _ < 

I COUNT 

WRITE (*, ' 

:A30, L3) 1 

) 

l 

*GETEDGE* 

STUCK = 1 , 

STUCK 

WRITE (* , ' 

( A3 0 , El 8 . 

6 ) 

• > ■ 

*GE TEDGE* 

YPREV = 1 , 

YPREV 

WRITE (*, ’ 

(ABO, 13) ' 

> 

i 

*GETEDGE * 

JU = ’ , 

JU 

WRITE (*, ' 

. A 3 0 , El 8 . 

6 ) 

1 ) ' 

* GETEDGE * 

Y ^ ‘ , 

Y 

WRITE (*, ' 

1 A3 C , 1 3 } ' 

) 

} 

*GETED GE* 

SSY - 1 , 

SSY 

WRITE (*, ' 

i A3 0 , L3 i 1 

) 

i 

* GETEDGE* 

SY.EQ.SSY - 

(SY.EQ.SSY) 

WRITE (*, ' 

( A3 0 , 613) 

’) 

i 

* GETEDGE* 

U2 :::: ' # 

(U2 (I) , 1 = 1, M) 

WRITE ( * , ' 

(ABC, 13) ’ 

) 

\ 

* GETEDGE * 

/ 

IX 

ENDIF 
IF (STUCK) 
PAUSE ( ’ S 

THEN 

Click in G 

~ i 


n Alt. DOWi 




ENDIF 

ENDDO 

ENDIF 

IF (SY.EQ.SSY) THEN 
JU = 0 
INFRONT = 0 . 

ISVERTEX = .FALSE. 

RETURN 

ENDIF 

XK = YPREV/ (YPREV-Y) 

CALL SETU (XU2 , U2 , U_MIN, U_MAX, M) 

Z2 = 0 . 

DO 1=1, M 

Z2 = Z2 +B1 (3, I) *XU2 (I) 

ENDDO 

Z1 = Z2 -U2 ( JU) *B1 ( 3 , JU) * (U MAX(JU)-U MIN(JUJ) 
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1198 Z = XK* Z2 + { 1 . -XK) *Z1 

1199 INFRONT = 1. 

1200 IF (Z.LT.O.) INFRONT = -1. 

1201 C 

1202 ISVERTEX = ( { Y . EQ . 0 . ) . AND . (ABS ( Z) . LT . TOL) ) 

1203 C 

1204 IF (DIAGS) THEN 

1205 WRITE ( * , ' : / A 5 C '■ • ) 1 Leaving GET EDGE ‘ 

1206 WRITE ( * , ' ; A3 0 , 6 I 3 :■ ' ) ' *GETSDGZ* U2 = ! , (U2 (I) , 1=1, M) 

1207 WRITE ( * , ' ; A 3 0 , 6 F 1 4 . 6 ) ’ ) 1 *GETEDGE* XU2 - ' , (XU2 (I) , 1 = 1, M) 

1208 WRITE ( * , * ■ A3 C , 1 3 } 1 ) ! *GSTEDGE* JU = ' , JU 

1209 WRITE (*, ' ' A3C, E13 . 6) ' ) ' * GET EDGE* XK = XK 

1210 WRITE (*, ' ;A30, E13 .6) ' ) 1 * GET EDGE* Z2 = Z2 

1211 WRITE (*, ' (A30, E18 . 6 i ' ) ' *GST2DGE* Z1 = ’ , Z1 

1212 WRITE ( * , ’ ' A 3 C , 213 .6 ■ ' ) ' * GET EDGE* Z - \ Z 

1213 WRITE (*, » :A30 # FIG .6; 1 ) ’ * GET EDGE* IN FRONT = ’ , INFRONT 

1214 WRITE ( * , ' A 3 : , 13 • 1 ) 1 *GF7EDGZ* LS7ER IEX - ' , ISVERTEX 

1215 END IF 

1216 C 

1217 RETURN 

1218 END 

1219 

1220 C 

1221 

1222 SUBROUTINE ISFACET ( ISOK, IUOUT, JUOUT, U123, 

1223 & IU, JU, Ul, U2 , B, M, TOL) 

1224 C 

1225 IMPLICIT NONE 

1226 

1227 LOG I CAL* 4 ISOK 

1228 INTEGER*4 IU, JU, M, Ul(20), U2{20), IUOUT, JUOUT, U123(20,3) 

1229 REAL *4 B(3,20), TOL 

1230 C LOCALS 

1231 INTEGER*4 UX1(20), UX2(20), I, J, K, II, JJ, KK 

1232 REAL*4 THEMAT (2,2) , MATINV{2,2), MATDET, T2(2), T1 (3 ) , TESTFACET ( 20 ) 

1233 REAL *4 TEMP 

1234 INTEGER*4 DIM, OBJ(20), UDEF(20), ITF(20), THEOBJ{20) 

1235 

1236 ISOK = .FALSE. 

1237 DO 1 = 1, M 

1238 UX1 (I) = Ul (I) 

1239 UX2 (I) = U2 (I) 

1240 ENDDO 

1241 Ul(IU) = 0 

1242 U2 ( JU) = 0 

1243 DIM = 0 

1244 DO I = 1 , M 

1245 OBJ ( I ) = U1(I) 

1246 IF ((I.EQ.IU) . OR . ( I . EQ . JU) .OR. (U2 { I ) . NE . Ul ( I ) ) ) THEN 

1247 OBJ ( I ) = 0 

1248 DIM = DIM + 1 

1249 UDEF (DIM) = I 

1250 END IF 

1251 ENDDO 

1252 C 

1253 IF (DIM . EQ . 2 ) THEN 

1254 ISOK = .TRUE. 
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1255 

1256 

1257 

1258 

1259 

1260 
1261 
1262 

1263 

1264 

1265 

1266 

1267 

1268 

1269 

1270 

1271 

1272 

1273 

1274 

1275 

1276 

1277 

1278 

1279 

1280 
1281 
1282 

1283 

1284 

1285 

1286 

1287 

1288 C 

1289 

1290 C 

1291 

1292 

1293 

1294 

1295 

1296 

1297 

1298 

1299 C 

1300 

1301 

1302 

1303 

1304 

1305 

1306 

1307 

1308 

1309 

1310 

1311 


JJ = MOD (II, 3) +1 
KK = MOD ( JJ , 3 ) +1 
THEMAT (1,1) = B ( I I , UDEF ( 1 ) ) 

THEMAT (1,2) = B { JJ, UDEF ( 1 ) ) 

THEMAT (2,1) = B ( II , UDEF (2) ) 

THEMAT (2,2) = 3 ( JJ, UDEF (2 ) ) 

MATDET = THEMAT (1,1)* THEMAT (2,2) - THEMAT (1,2) * THEMAT (2,1) 

IF (MATDET. NE. 0 . ) THEN 

MATINV (1,1) = THEMAT (2, 2) /MATDET 
MAT I NV (1,2) = - THEMAT (1, 2) /MATDET 
MATINV (2, 2) = THEMAT (1, 1) /MATDET 
MATINV (2,1) = - THEMAT ( 2, 1) /MATDET 

T2 (1) = -MATINV ( 1 , 1) *B (KK, UDEF (1) ) -MATINV (1, 2) *B (KK, UDEF (2 ) ) 

T2 (2) = -MATINV (2 , 1) *B ( KK, UDEF ( 1 ) ) -MATINV (2 , 2) *B (KK, UDEF ( 2 ) ) 
Tl(KK) = 1. 

T1 (II) = T2 (1) 

T1 (JJ) = T2 (2) 

DO 1 = 1 , M 

TESTFACET ( I ) = 0. 

ITF(I) = 0 
UDEF ( I ) = 0 
DO J=1 , 3 

TESTFACET ( I ) = TESTFACET ( I ) +T1 (J) *B (J , I ) 

ENDDO ! DO J=l,3 

IF (ABS (TESTFACET (I) ) .LT.TOL) THEN 
ITF(I) = 0 

ELSEIF (TESTFACET ( I) . GT . 0 . ) THEN 
ITF(I) = 1 

ELSEIF (TESTFACET (I) . LT . 0 . ) THEN 
ITF(I) = -1 
ENDIF 

ENDDO ! DO 1 = 1, M 

DIM = 0 
DO 1=1, M 

THEOBJ(I) = OBJ (I) 

IF ( (OBJ ( I ) . EQ . 0 ) .OR. (ITF(I) .EQ.0) .OR. ( ITF ( I ) . NE . OB J ( I ) ) ) 
THEOBJ(I) = 0 
DIM = DIM + 1 
UDEF (DIM) = I 
ENDIF 

ENDDO ! DO 1 = 1, M 

IF ( DIM . NE . 2 ) THEN 
DIM = 0 
J = 0 
DO 1=1, M 

THEOBJ ( I ) = OBJ (I) 

IF ( (OBJ(I) .EQ.0) .OR. ( ITF { I ) . EQ . 0 ) .OR. ( - ITF ( I ) . NE . OB J ( I ) 
THEOBJ (I) = 0 
DIM = DIM + 1 
UDEF (DIM) = I 
ENDIF 

ENDDO I DO 1 = 1, M 
ENDIF 1 IF (DIM . NE . 2 ) THEN 


THEN 


) ) THEN 
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1312 C 


1313 

IF (DIM.NE.2) ISOK = .FALSE. 

1314 


1315 

GOTO 194 

1316 


1317 

END IF ! IF (MATDET.NE.0.) THEN 

1318 

ENDDO ! DO 11=1,3 

1319 194 

CONTINUE 

1320 

END IF i IF (DIM.2Q.2) 

1321 C 


1322 

IF (ISOK) THEN 

1323 

IUOUT = UDEF(l) 

1324 

JUOUT = UDEF ( 2 ) 

1325 

DO 1=1, M 

1326 

DO J=l, 3 

1327 

U12 3 (I, J) = OBJ ( I ) 

1328 

ENDDO 

1329 

ENDDO 

1330 

U123 (IUOUT, 1) = 1 

1331 

U123 (JUOUT, 1) = 1 

1332 

U123 (IUOUT, 2) = -1 

1333 

U123 (JUOUT, 2) = 1 

1334 

U123 (IUOUT, 3) = 1 

1335 

U123 (JUOUT, 3) = -1 

1336 

ELSE 

1337 

IUOUT = 0 

1338 

JUOUT = 0 

1339 

DO 1=1, M 

1340 

DO J=l, 3 

1341 

U12 3 (I, J) = 0 

1342 

ENDDO 

1343 

ENDDO 

1344 

END IF 

1345 


1346 

RETURN 

1347 

END 

1348 


1349 C 


1350 


1351 

SUBROUTINE SORTC (X, IY, N, XS , IYC) 

1352 


1353 

IMPLICIT NONE 

1354 


1355 

INTEGERM I, IL_SORTC ( 36 ) , IMED, IP1 

1356 

INTEGERM IY(20), IYC(20), J, JMI, Ji 

1357 

INTEGER* 4 NM1 

1358 

REAL* 4 HOLD, AMED, TX, X(20), XS(20) 

1359 C 


1360 C 

CHECK THE INPUT ARGUMENTS FOR ERRORS 

1361 C 


1362 

IPR=11 

1363 

I F ( N . LT . 1 ) GOTO 5 0 

1364 

IF (N . EQ . 1 ) GOTO 5 5 

1365 

HOLD=X ( 1 ) 

1366 

DO60I=2 , N 

1367 

IF (X (I ) . NE . HOLD ) GOTO 90 

1368 60 

CONTINUE 


ITY, IU_S0RTC (36) 
L, LMI , M, MID, N 
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1369 C WRITEt*, 9 ) HOLD 

1370 D06 1 I = 1 , N 

1371 XS(I)=X{I) 

1372 I YC ( I ) = IY ( I ) 

1373 61 CONTINUE 

1374 RETURN 

1375 C 50 WRITE (*,15) 

1376 C WRITE ( * , 4 7 ) N 

1377 50 RETURN 

1378 C 55 WRITE (* , 18 ) 

1379 55 XS(1)=X(1) 

1380 IYC (1) =IY (1) 

1381 RETURN 

1382 90 CONTINUE 

1383 9 FORMAT (IX ,'***** NON- RATAL DIAGNOS 

1384 1NT A VECTOR) TO THE SCETC SUBROUTINE HAS ALL ELEMENTS = ' , E15 . 3 , 

Q *T. 1 '*'*■*★* i J 

1386 15 FORMAT (IX , '***«* FATAL ERROR- -THE SECOND INPUT ARGUMENT TO THE 

13 3 7 : 50 RTC SUBROUTINE IS NON-POSITIVE ***»» ' ) 

1388 18 FORMAT ( IX , ’****♦ NON -FATAL DIAGNOSTIC- -THE SECOND INPUT ARGUME 

13 89 l • J i Hi ioORTL SU3Roo i INi HAS i'HE VaLOii 1 * * * * * 1 ) 

1390 47 FORMAT (IX , ****** THE VALUE OF THE ARGUMENT IS * , 18 f ’ ***** 1 ) 

1391 C 

1392 C START POINT 

1393 C 

1394 C COPY THE VECTOR X INTO THE VECTOR XS 

1395 DO100I=l,N 

1396 XS(I)-X(I) 

1397 100 CONTINUE 

1398 C 

1399 C COPY THE VECTOR IY INTO THE VECTOR IYC 

1400 C 

1401 DO150I«l,N 

1402 IYC(I) =IY (I) 

1403 150 CONTINUE 

1404 C 

1405 C CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED 

1406 C 

1407 NM1=N- 1 

1408 DO200I=l , NM1 

1409 IP1=I+1 

1410 IF(XS (I) .LE.XS (IPX) ) GO TO200 

1411 GOTO250 

1412 200 CONTINUE 

1413 RETURN 

1414 250 M«1 

1415 1=1 

1416 J=N 

1417 305 I F ( I . GE . J ) GOT03 7 0 

1418 310 K= I 

1419 MID- ( I + J) /2 

1420 AMED=XS (MID) 

1421 I MED = IYC (MID) 

1422 IF(XS(I) .LE. AMED) GOTO320 

1423 XS (MID) =XS (I) 

1424 IYC (MID) = I YC ( I ) 

1425 XS ( I ) =AMED 
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1426 

1427 

1428 

1429 

1430 

1431 

1432 

1433 

1434 

1435 

1436 

1437 

1438 

1439 

1440 

1441 

1442 

1443 

1444 

1445 

1446 

1447 

1448 

1449 

1450 

1451 

1452 

1453 

1454 

1455 

1456 

1457 

1458 

1459 

1460 

1461 

1462 

1463 

1464 

1465 

1466 

1467 

1468 

1469 

1470 

1471 

1472 

1473 

1474 

1475 

1476 

1477 

1478 

1479 

1480 

1481 

1482 


IYC(I) =IMED 
AMED=XS (MID) 

IMED= I YC (MID) 

320 L= J 

I F (XS ( J) . GE . AMED ) G0T03 4 0 
XS (MID) = XS ( J) 

I YC (MID) = I YC (J) 

XS ( J) =AMED 
I YC (J) =IMED 
AMED-XS (MID) 

IMED=IYC (MID) 

IF (XS ( I ) . LE . AMED ) G0T03 4 0 
XS (MID) =XS (I) 

I YC (MID) =IYC ( I ) 

XS (I) =AMED 
IYC(I) =IMED 
AMED=XS (MID) 

IMED=I YC (MID) 

GOTO340 

330 XS (L) =XS (K) 

I YC (L) =1 YC ( K) 

XS (K) =TX 
IYC(K) =ITY 

340 L=L- 1 

IF (XS (L) .GT.AMED) GOTO340 
TX=XS (L) 

ITY = I YC ( L) 

350 K=K+1 

IF (XS { K) . LT . AMED ) G0T03 5 0 
IF (K. LE . L) GOTO330 
LMI=L- I 
JMK= J-K 

IF (LMI .LE. JMK)GOTO360 
IL_SORTC(M) =1 
IU_SORTC(M) =L 
I =K 
M=M+ 1 
GOTO 3 8 0 

360 IL_SORTC (M) =K 
IU_SORTC (M) =J 
J=L 
M = M+1 
GOTO380 

370 M=M - 1 

IF (M.EQ. 0) RETURN 
I=IL_SORTC(M) 

J=IU_SORTC (M) 

380 JMI=J - I 

IF ( JMI . GE . 11) GOTO310 
IF (I . EQ.l)GOTO305 
1 = 1-1 

390 1=1+1 

IF ( I . EQ . J) GOTO370 
AMED=XS ( 1+1) 

IMED=IYC ( 1+1) 

I F (XS ( I ) . LE . AMED) GOTO390 
K= I 
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1483 395 XS (K+l) =XS (K) 

1484 I YC (K+ 1 ) = I YC ( K) 

1485 K=K- 1 

1486 I F ( AMED . LT . XS (K) )G0T0395 

1487 XS ( K+ 1 ) =AMED 

1488 I YC (K+l) =IMED 

1489 GOTO 390 

1490 

1491 RETURN 

1492 END 

1493 

1494 C 

1495 

1496 SUBROUTINE M I NNORM (UM INNORM , SCALE, 

1497 C AS A FUNCTION OF 

1498 & P, U, UTRIM, MTRIM, UDA, B, U_MIN, U_MAX, M, TIME, XSCALE) 

1499 

1500 IMPLICIT NONE 

1501 C ** Parameters 

1502 

1503 

1504 C ** INPUTS: 

1505 

1506 REALM P(20,3), U (20) , UTRIM(20), UDA(20) 

1507 REALM B(3,20), MTRIM ( 3 ) , U_MAX{20), U_MIN(20) 

1508 REALM TIME, XSCALE 

1509 INTEGERM M, IMODE 

1510 

1511 C ** OUTPUTS: 

1512 

1513 REALM UMINNORM ( 2 0 ) , SCALE 

1514 

1515 C ** LOCALS: 

1516 

1517 REALM UKDA(20), UP(20), UDELTA ( 20 ) , M_ATT ( 3 ) 

1518 INTEGERM I, J, K 

1519 C 

1520 SCALE = 1. 

1521 

1522 DO 1=1, M 

1523 UKDA(I) = U(I) + UDA(I) 

1524 ENDDO 

1525 

1526 DO 1 = 1,3 

1527 M_ATT(I)=0. 

1528 DO J=1 , M 

1529 M_ATT ( I ) =M_ATT(I) +B ( I , J) *UKDA(J) 

1530 ENDDO 

1531 ENDDO 

1532 

1533 DO 1 = 1, M 

1534 IF ( (U_MIN(I) .EQ.0. ) .OR. (U_MAX(I) . EQ . 0 . ) ) SCALE = 0. 

1535 UKDA(I) = U(I) + UDA(I) 

1536 UP(I) = UTRIM (I) 

1537 DO J=l, 3 

1538 UP(I) = UP(I) + P ( I , J) * (M_ATT ( J) -MTRIM ( J) ) 

1539 ENDDO 
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1540 

1541 

1542 

1543 

1544 

1545 

1546 

1547 

1548 

1549 C 

1550 

1551 C 

1552 

1553 

1554 

1555 C 

1556 

1557 

1558 

1559 C- 


UDELTA ( I ) = SCALE* (UP (I) -UKDA(I) ) 

IF ( UDELTA ( I ) . NE . 0 . ) THEN 

IF ( (UDELTA (I) . GT . U_MAX ( I ) ) .AND. (U_MAX(I) .GT.O. ) ) 

SCALE = U_MAX ( I ) /UDELTA ( I ) 

ELSEIF ( (UDELTA ( I ) .LT.U_MIN(I) ) .AND. (U_MIN(I) . LT . 

SCALE = U_MIN ( I ) /UDELTA ( I ) 

END IF 
END IF 
ENDDO 

SCALE = XSCALE* SCALE 
DO 1=1, M 

UMINNORM ( I ) = UKDA(I) + SCALE* (UP ( I ) -UKDA ( I ) ) 

ENDDO 

RETURN 

END 


THEN 

0 . ) ) THEN 



